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The  major  objective  of  this  study  is  to  analyze  buckling  and  delamination 
behavior  of  composite  stiffened  panels  subjected  to  axial  compression. 

First,  a  combined  analytical  and  experimental  study  of  a  blade  stiffened  composite 
panel  subjected  to  axial  compression  was  conducted.  The  effects  of  the  differences 
between  a  simple  model  used  to  design  the  panel  and  the  actual  experimental  conditions 
were  examined.  It  was  found  that  in  spite  of  many  simplifying  assumptions  the  design 
model  did  reasonably  well  in  that  the  experimental  failure  load  was  only  10%  higher  than 
the  design  load.  Several  structural  analysis  programs,  including  PANDA2,  STAGS,  and 
ABAQUS,  were  used  to  obtain  high  fidelity  analysis  results.  The  buckling  loads  from 
STAGS  agreed  well  with  the  experimental  failure  loads.  However,  substantial  differences 
were  found  in  the  out-of-plane  displacements  of  the  panel.  Efforts  were  made  to  identify 
the  source  of  these  differences.  Implementing  non-uniform  load  introduction  with  general 
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contact  definition  in  the  STAGS  finite  element  model  improved  correlation  between  the 
measured  and  predicted  out-of-plane  deformations. 

Next,  a  new  method  called  Crack  Tip  Force  Method  (CTFM)  is  derived  for 
computing  point-wise  energy  release  rate  along  the  delamination  front  in  delaminated 
plates.  The  CTFM  is  computationally  simple  as  the  G  is  computed  using  the  forces 
transmitted  at  the  crack-tip  between  the  top  and  bottom  sub-laminates  and  the  sub- 
laminate  properties. 

Finally,  buckling  and  postbuckling  of  a  blade-stiffened  composite  panel  under 
axial  compression  with  a  partial  skin-stiffener  debond  are  investigated.  Two  different 
finite  element  models,  where  nodes  of  the  panel  skin  and  the  stiffener  flange  are  located 
on  the  mid-plane  or  at  the  interface  between  skin  and  flange,  are  used.  Linear  buckling 
analysis  is  conducted  using  both  STAGS  and  ABAQUS.  Postbuckling  analysis  is 
conducted  with  STAGS.  Comparison  between  the  present  results  and  previous  buckling 
analysis  results  show  a  good  correlation.  Buckling  analysis  results  for  various  stiffener 
geometries  and  debond  ratios  are  presented. 
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CHAPTER  1 
INTRODUCTION 

1.1  Background 

Stiffened  laminated  composite  panels  have  been  considered  for  use  in  weight- 
sensitive  structures  such  as  aircraft  and  missile  structural  components.  The  main 
advantage  of  the  stiffeners  is  the  increased  structural  efficiency  of  the  structure  with  a 
minimum  of  additional  material.  Due  to  the  high  stiffness  of  fiber  composites,  stiffened 
composite  panels  are  usually  thin.  Thus,  buckling  characteristics  are  critical 
considerations  for  the  optimum  design  of  composite  structures  made  of  laminated 
composite  plates.  Buckling  depends  on  a  variety  of  factors,  such  as  the  geometry  of  the 
members,  boundary  conditions  and  material  properties.  Because  of  geometric 
complexity,  stiffened  laminated  composite  panels  in  axial  compression  have  several 
buckling  modes  including  general  instability,  local  skin  buckling  and  rolling  of 
stiffeners.  Furthermore,  stiffened  laminated  composite  flat  panels  usually  exhibit  stable 
postbuckling  behavior  which,  in  general,  leads  to  significant  differences  between 
buckling  load  and  ultimate  failure  load.  The  correct  estimation  of  the  load  carrying 
capacity  of  stiffened  panels  is  therefore  very  complicated. 

Advanced  fiber-reinforced  composite  materials  such  as  graphite-epoxy  have 
relatively  low  transverse  tensile  and  interlaminar  shear  strengths  compared  to  in-plane 
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strength.  Therefore,  stiffened  panels  made  of  laminated  composite  are  very  susceptible 
to  delamination.  Delamination,  also  referred  to  as  interface  cracking  or  debonding  where 
adjacent  laminae  become  separated  from  one  another,  has  been  one  of  the  major  known 
weaknesses  of  laminated  composite  structures.  Delamination,  coming  from  initial 
manufacturing  imperfections  or  in-service  damage  such  as  foreign  object  impact,  can 
significantly  reduce  the  stiffness  and  strength  of  the  composite  structures.  It  is  also  well 
known  that  delamination  and  skin-stiffener  separation  are  common  failure  modes  of  a 
stiffened  composite  panel  in  axial  compression.  When  a  delamination  is  present,  it  is 
very  important  to  identify  not  only  its  global  influence  on  the  load  carrying  capacity  of 
the  structure  but  also  its  local  behavior  under  the  applied  load.  Energy  release  rate  has 
been  accepted  as  a  measure  for  predicting  delamination  propagation.  Most  available 
methods  of  computing  energy  release  rate  use  2-D  or  3-D  solid  finite  elements.  However 
it  is  computationally  very  expensive  to  use  solid  elements  for  modeling  entire 
complicated  structures  made  of  laminated  composite  materials.  Thus  simplified  beam, 
plate  and  shell  theories  are  frequently  used  for  structural  analysis  of  those  structures. 
Therefore,  it  may  be  a  good  idea  to  also  use  the  same  theory  to  calculate  energy  release 
rate. 

1.2  Literature  Survey 
1-2.1  Buckling  and  Postbuckling  Analysis  of  Stiffened  Panels 

Research  on  the  buckling  and  postbuckling  behavior  of  stiffened  panels  has  been 
of  interest  for  many  years,  with  many  researchers  exploring  the  response  of  the  stiffened 
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panels.  However,  due  to  their  geometric  complexity  and  the  many  parameters  involved, 
a  complete  understanding  of  all  aspects  of  behavior  is  not  yet  fully  achieved.  Several 
researchers  [1-5]  compiled  surveys  on  buckling  and  postbuckling  behavior  of  composite 
panels.  Leissa  [1]  compiled  extensive  results  on  buckling  and  postbuckling  behavior  of 
laminated  composite  panels.  Analytical  effects  of  various  boundary  conditions,  stacking 
sequences,  and  transverse  shear  deformation  on  buckling  and  postbuckling  behavior 
were  explored.  Noor  and  Peters  [2]  reviewed  two  aspects  of  the  numerical  simulation  of 
the  buckling  and  postbuckling  responses  of  composites  structures.  The  first  aspect  was 
exploiting  non-traditional  symmetries  exhibited  by  composite  structures,  and  strategies 
for  reducing  the  size  of  the  model  and  the  cost  of  buckling  and  postbuckling  analyses  in 
the  presence  of  symmetry-breaking  conditions  (e.g.,  asymmetry  of  the  material, 
geometry,  and  loading).  The  second  aspect  was  the  prediction  of  onset  of  local 
delamination  in  the  postbuckling  range  and  accurate  determination  of  transverse  shear 
stresses  in  the  structure.  Bushnell  [3]  divided  the  literature  in  the  field  of  buckling  of 
stiffened  panel  into  three  categories.  One  in  which  structural  analysis  is  emphasized,  a 
second  in  which  optimum  design  is  emphasized,  and  a  third  in  which  test  results  are 
emphasized.  Recently,  Knight  and  Stames  [4]  reviewed  some  of  the  historic 
developments  of  shell  buckling  analysis  and  design  and  identified  key  research 
directions  for  reliable  and  robust  methods  in  shell  stability  analysis  and  design.  Bedair 
[5]  presented  an  extensive  literature  review  on  stability  of  stiffened  panel  under  uniform 
compression.  He  classified  the  literature  into  two  categories,  analysis  and  design.  The 
objective  of  the  first  category  is  to  develop  numerical  or  analytical  formulations  to 
predict  the  global  and  local  buckling  load  of  structures.  In  doing  so,  several  assumptions 


are  postulated  in  idealizing  the  structure  in  order  to  facilitate  a  solution.  The  objective  of 
the  second  group  of  researchers  is  to  develop  simplified  models  to  predict  the  ultimate 
strength,  or  collapse  load  of  the  structures.  Several  simplified  models  have  been 
developed  for  that  purpose. 

The  main  objective  of  this  section  is  to  review  briefly  previous  works  in  the  area 
of  buckling  and  postbuckling  behavior  of  stiffened  composite  panels.  Surveys  are 
divided  into  two  categories,  analysis  methods  and  correlation  of  analysis  and 
experiment. 

1.2.1.1  Analysis  methods 

The  analysis  of  stiffened  panel  can  be  performed  by  either  smeared  or  discrete 
stiffener  approach.  A  smeared  stiffener  approach  converts  the  stiffened  plate  into  an 
equivalent  plate  with  constant  thickness  by  smearing  out  the  stiffeners.  This  method 
provides  accurate  analysis  results  of  global  buckling  load  of  stiffened  panels  when 
stiffeners  are  identical  and  closely  spaced.  However,  this  approach  ignores  the  discrete 
nature  of  the  structures  and  does  not  consider  all  potential  buckling  modes.  The  discrete 
stiffener  approach  does  not  have  limitation  of  stiffener  spacing  and  uniformity. 

In  early  studies  of  buckling  analysis  of  stiffened  panels,  Wittrick  and  William  [6] 
developed  a  general -purpose  computer  program,  VIP  ASA,  for  determining  the  critical 
buckling  stresses  and  natural  frequencies  of  thin  prismatic  structures  that  consist  of  a 
series  of  flat  plate  connected  rigidly  along  their  longitudinal  edges.  The  response  of  each 
plate  element  making  up  the  stiffened  panel  is  obtained  from  an  exact  solution  of  thin- 
plate  theory.  This  approach  is  commonly  referred  to  as  exact  finite  strip  method.  The 
analysis  used  in  VIP  ASA  is  similar  to  that  which  Viswanathan  et  al.  [7]  incorporated  in 


their  computer  program  BUCLASP2.  BUCLASP2  assumes  that  panel  elements  are 
orthotropic  and  have  balanced  laminates,  material  is  linear  elastic  and  thin-plate  theory. 
Stroud  and  Agranoff  [8-9]  proposed  a  simplified  analysis  based  on  buckling  of 
othotropic  plates  with  simply  supported  boundary  conditions.  The  global  buckling 
analysis  of  stiffened  panel  was  conducted  as  an  orthotropic  plate  with  smeared 
stiffeners,  assuming  as  a  wide  column.  In  contrast  with  the  simplified  analysis,  the 
VIPASA  analysis  code  provides  a  high-quality  buckling  analysis  that  considers  all 
buckling  modes  and  ensures  continuity  of  the  buckling  pattern  across  the  intersection  of 
neighboring  plate  elements.  Stroud  and  Anderson  [10]  developed  a  stiffened  panel 
design  code  PASCO,  which  uses  VIPASA  analysis  as  well  as  design  techniques  from 
the  early  simplified  methods  and  expanded  capabilities  such  as  initial  bow-type 
imperfection,  bending  moments,  and  temperature  loading.  The  major  known  drawback 
of  the  VIPASA  code  is  underestimation  of  the  shear  buckling  modes  when  a  buckling 
half-wave  length  is  equal  to  the  panel  length  [10,12].  To  overcome  inaccuracy  involving 
shear  load  and  anisotropy  exhibited  by  VIPASA  and  PASCO,  Anderson  et  al.  [11] 
developed  the  computer  program  VICON.  The  analysis  of  VICON  assumes  that  the 
deflection  of  the  plate  assembly  can  be  expressed  as  a  Fourier  series,  which  can  be  used 
to  calculate  the  forces  at  the  longitudinal  junctions  between  the  plates  by  the  same 
stiffness  matrices  that  result  from  the  VIPASA  analysis.  The  total  energy  of  the  panel  is 
expressed  in  terms  of  VIPASA  stiffness  matrices  plus  conventional  stiffness  of  the 
supporting  structure.  Then  the  total  energy  is  minimized  subject  to  the  constraints,  using 
Lagrange  multipliers,  to  obtain  the  governing  equations. 

Bushnell  [13-14]  developed  PANDA,  an  early  version  of  PANDA2,  where 


buckling  loads  are  computed  by  the  use  of  simple  assumed  displacement  functions  used 
in  conjunction  with  Donnell-type  kinematic  relation.  Several  types  of  general  and  local 
buckling  modes  were  included.  VIPASA,  PASCO,  and  PANDA  cannot  perform 
nonlinear  postbuckling  analysis  of  the  stiffened  panel.  Bushnell  [15-17]  released 
PANDA2,  which  incorporated  nonlinear  theory  for  prediction  of  behavior  of  locally 
imperfect  panels.  Nonlinear  strain-displacement  relations  analogous  to  those  developed 
in  1946  by  Koiter  for  perfect  panels  were  extended  to  handle  panels  with  imperfections 
in  the  form  of  critical  local  bifurcation  buckling  mode.  In  PANDA2,  local  and  general 
buckling  loads  are  calculated  with  use  of  either  closed-form  expression  [14]  or  with  use 
of  discretized  models  of  panel  cross-sections  [15-17].  The  discretized  model  is  based  on 
one-dimensional  discretization  that  is  commonly  referred  to  as  approximate  finite  strip 
method.  Approximate  finite  strip  method  assumes  displacement  of  strip  can  be 
expressed  with  combination  of  crosswise  polynomial  function  and  lengthwise 
trigonometric  function.  Dowe  and  his  group  [18-22]  have  used  the  method  in 
conjunction  with  a  non-linear  analysis  based  on  first-order  shear  deformation  theory. 
The  Finite  strip  method,  which  lies  between  the  conventional  Rayleigh-Ritz  method  and 
the  finite  element  method,  provide  a  means  of  solving  prismatic  plate-structure  problems 
with  an  attractive  blend  of  accuracy,  economy  and  ease  of  modeling  [22]. 

The  finite  element  method  is  the  most  powerful  method  to  predict  the  buckling 
and  postbuckling  behavior  of  structures.  The  historic  developments  of  shell  buckling 
analysis  using  finite  element  up  to  1997  can  be  found  in  Knight  and  Stames  [5].  With 
advancing  computer  power,  the  finite  element  method  is  widely  used  to  solve  various 
engineering  problems.  Many  current  commercially  available  finite  element  codes  such 
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as  MSC/NASTRAN  [23]  and  HKS/ABAQUS  [24]  provide  buckling  and  postbuckling 
analysis  capability  of  the  stiffened  composite  panels  as  one  of  several  analysis  options. 
Knight  and  Stames  [5]  pointed  out  in  their  review  paper  that  STAGS  [25]  is  perhaps  the 
premier  shell  analysis  code  that  focused  primarily  on  shell  analysis  and  solution 
procedures  for  shell  problems. 
1.2.1.2  Analytical  and  experimental  correlation 

Williams  and  Stein  [26]  examined  J-  and  blade-stiffened  graphite/epoxy  panels 
experimentally  as  well  as  analytically  using  several  analysis  codes  such  as  VIPASA, 
BUCLASP-2  and  an  early  version  of  STAGS,  which  was  based  on  a  two  dimensional 
finite  difference  technique.  In  their  study,  correlation  of  experimental  and  analytical 
results,  which  included  inplane  displacement  restraints,  indicated  that  the  buckling  strain 
of  J-stiffened  specimens  were  75%  to  80%  of  analytical  values  and  that  of  blade- 
stiffened  panels  were  84%  to  97%  of  the  analytical  values.  A  nonlinear  response  was 
exhibited  by  several  of  the  specimens  in  which  large  lateral  displacement  in  the  order  of 
one-quarter  of  the  thickness  of  the  panel  plate  segments  were  observed. 

Stames  and  coworkers  [27-28]  investigated  the  postbuckling  behavior  of  flat  and 
curved  stiffened  graphite-epoxy  panels  loaded  in  compression.  Panels  with  four  equally 
spaced  I-shaped  stiffeners  and  quasi-isotropic  skin  were  tested.  Failure  of  all  panels 
initiated  in  a  skin-stiffener  interface  region.  They  showed  that  analytical  results  obtained 
from  PASCO  as  well  as  STAGS  correlate  well  with  typical  postbuckling  test  results  up 
to  failure.  Their  results  also  showed  that  modeling  of  the  stiffener  components  with  plate 
elements  having  appropriate  stiffness  is  required  to  obtain  satisfactory  correlation  with 
the  postbuckling  test  results. 
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Romeo  [29]  conducted  several  tests  on  graphite-epoxy  hat-and  blade-stiffened 
panels  under  uniaxial  compression  and  wing-box  beams  under  pure  bending  to  verify  the 
accuracy  of  the  theoretical  analysis.  Overall  buckling,  local  buckling  and  torsional 
buckling  were  determined  separately  using  a  simple  engineering  formula,  and 
interactions  between  these  modes  were  not  considered.  Adequate  correlation  with 
experimental  results  was  obtained  for  axial  compression  when  the  Euler  or  torsional 
buckling  mode  was  critical;  buckling  occurred  at  lower  strain  values  than  predicted 
when  the  local  buckling  mode  was  critical.  Furthermore,  he  showed  that  simple 
compression  tests  could  not  represent  the  load  conditions  of  wing-box  compression 
panel  properly;  in  particular,  the  bending  curvature  causes  a  distributed  load 
perpendicular  to  panels  that  could  reduce  the  longitudinal  load  at  which  buckling 
occurred. 

Bushnell  et  al.  [30]  conducted  optimum  design,  fabrication,  and  test  of  graphite- 
epoxy  curved,  locally  buckled  panels  in  axial  compression.  Three  nominally  identical 
large  panels  were  tested.  Two  of  three  tests  gave  reasonably  good  agreement  between 
test  and  theory,  both  with  regard  to  loads  at  which  the  panels  failed  and  the  mode  of 
failure.  They  also  conducted  experimental  comparison  between  specimens  with  stitched 
skin-flange  combination  and  specimens  with  adhesive-bonded  skin-flange  combination. 
They  found  that  load  carrying  capacity  of  stitched  specimens  were  lower  than  those  of 
adhesively  bonded  specimens. 

Wieland  et  al.  [31]  investigated  the  buckling,  postbuckling  and  crippling  of 
AS4/3502  graphite-epoxy  Z-section  stiffeners  as  a  function  of  specimen  structural 
parameters.  Variables  considered  were  flange  and  web  widths,  flange-to-web  comer 


parameters.  Variables  considered  were  flange  and  web  widths,  flange-to-web  comer 
radius,  thickness  and  stacking  sequence.  Analytical  model  was  based  on  a  classical 
model  on  local  buckling  and  the  numerical  analyses  were  conducted  with  the  ABAQUS 
finite  element  code.  The  results  showed  that  the  nature  of  the  load  redistribution  after 
buckling  and  its  effect  on  postbuckling  stiffness  are  related  to  the  geometric  variables. 
The  agreement  between  the  analytical  and  experimental  buckling  loads  was  generally 
good.  The  agreement  degrades  as  the  ratio  of  flange  width  with  respect  to  web  width  is 
reduced  and  as  the  section  comer  radius  is  increased. 

Fan  et  al.  [32]  performed  the  pre-  and  post-buckling  analysis  for  stiffened  panels. 
Both  the  thick-wall  stiffeners  with  square  cross  section  and  the  thin-wall  blade  stiffeners 
were  employed  in  their  study.  After  linear  buckling  analysis,  they  used  an  incremental 
analysis  along  the  load  path  with  a  special  iteration  technique,  called  initial  value 
method  to  improve  the  normal  Newton-Raphson  method.  The  computational  results  of 
both  displacement  control  and  load  control  were  presented.  The  computed  local  buckling 
load  and  buckling  mode  agreed  well  with  the  test  results.  However,  the  computed  global 
buckling  load  with  uniform  end-shortening  displacement  was  much  higher  than  the  test 
results.  But,  the  computed  global  buckling  load  with  load  control  was  close  to  the  test 
results.  The  reason  of  this  discrepancy  was  not  clearly  explained.  Three  failure  criteria  of 
composite  stiffened  panel,  which  are  maximum  strain,  debond  failure,  and  combined 
global-local  buckling  criterion  were  proposed. 

Nagendra  et  al.  [33]  studied  the  optimum  design  of  blade  stiffened  panels  with 
holes  under  buckling  and  strain  constraints.  They  used  PASCO  for  buckling  analysis  and 
optimization  with  continuous  thickness  design  variables  and  the  Engineering  Analysis 
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Language  (EAL)  finite  element  analysis  code  [34]  for  calculating  strain  and  their 
derivatives  with  respect  to  design  variables.  Later  their  optimally  designed  panels  with 
and  without  centrally  located  holes  were  tested  and  analytical  and  experimental  results 
were  compared  focusing  prebuckling  behavior  of  the  panel  [35].  Prebuckling  stiffness 
from  test  were  about  10%  lower  than  analytical  values  and  failure  loads  from  test  was 
also  10%  lower  than  that  from  design.  Since  the  uncertainties  in  the  geometric  and 
material  properties  did  not  account  for  the  discrepancy  between  analytical  and 
experimental  buckling  loads,  they  hypothesized  that  geometric  imperfections  and 
eccentricities  may  had  reduced  the  buckling  load. 

Chow  and  Atluri  [36]  proposed  failure  criterion  of  mixed-mode  stress  intensity 
factors  for  the  postbuckling  strength  of  stiffened  panel.  They  showed  that  post-buckling 
strength  of  the  stiffened  panels  compare  quite  favorably  with  the  experimental  results  of 
Stames  et  al.  [27]  and  the  standard  deviation  of  the  error  was  less  than  10  %. 

Young  and  Hyer  [37]  presented  modeling  procedures  that  predict  the 
postbuckling  response  of  composite  panels  with  skewed  stiffeners.  Five  panel 
configurations  with  various  combinations  of  skin  and  stiffener  orientation  were  tested.  A 
uniform  end  shortening  displacement  was  applied  to  the  upper  end  of  the  panel  in  the 
axial  direction,  and  the  axial  displacement  of  the  lower  end  was  restrained.  The  upper 
and  lower  ends  were  clamped  and  the  unloaded  sides  were  simply  supported.  The 
individual  effect  of  shell  elements,  potted  load  introduction,  material  properties,  and 
initial  geometric  imperfections  was  examined.  The  results  showed  that  inaccurate 
modeling  assumptions  and  anomalies  in  the  test  such  as  the  support  fixtures,  the  loading 
frame,  and  the  load  introduction  of  the  test  specimen  could  cause  the  predicted  response 
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and  the  measured  response  to  differ  substantially. 

Davila  et  al.  [38]  conducted  progressive  failure  analysis  for  the  simulation  of 
damage  initiation  and  growth  in  stiffened  thick-skin  stitched  graphite-epoxy  panels 
loaded  in  axial  compression.  Failure  indices  approach,  proposed  by  Chang  and  Chang 
[39],  was  adopted  to  evaluate  the  failure  mode  and  location  corresponding  to  all  of  the 
major  composite  laminate  failure  modes  except  delamination.  Superposed  layers  of  shell 
elements  with  multiple  integration  points  through  the  thickness  were  used  to  separate  the 
failure  modes  for  each  ply  orientation  and  to  obtain  the  correct  effect  of  bending  loads 
on  damage  progression.  The  analysis  results  were  compared  with  experimental  results 
for  three  stiffened  panels  with  notches  oriented  at  0,  15,  and  30  degrees  to  the  panel 
width  dimension  and  found  to  be  in  excellent  correlation  with  the  experimental  results. 
The  local  reinforcing  effect  of  Kevlar  stitches  was  simulated  in  the  finite  element  model 
by  multiplying  the  fiber  buckling  strength  allowable  value,  independent  of  the  other 
stress  components,  by  a  stitch  factor  that  is  determined  empirically.  A  parametric  study 
was  performed  to  investigate  the  damage  growth  retardation  characteristics  of  the  Kevlar 
stitch  lines  in  the  panels.  The  debond  between  the  stiffener  flange  and  the  skin  were  not 
modeled.  Hence,  the  predicted  results  were  found  to  be  less  accurate  after  the  damage 
zone  reached  the  stiffener  flange. 

Singer  et  al.  [40]  presented  conventional  and  less  conventional  experimental 
methods  in  buckling  of  a  vast  variety  of  thin-walled  structures  in  considerable  detail. 
The  parameters,  which  may  influence  the  test  results,  were  systematically  highlighted: 
imperfections,  boundary  conditions,  loading  conditions,  effect  of  holes  and  cutouts. 
Though  authors  deals  primarily  with  experimental  methods  and  test  results,  the 
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theoretical  concepts  of  the  basic  instabihty  phenomena  and  numerical  methods  were  also 
briefly  reviewed. 

Sleight  [41]  analyzed  a  composite  blade-stiffened  panel  with  a  discontinuous 
stiffener  loaded  in  axial  compression.  A  progressive  failure  analyses  using  Hashin's 
criterion  [42]  was  performed  on  the  blade-stiffened  panel.  The  progressive  failure 
analysis  and  test  results  showed  good  correlation  up  to  the  load  where  local  failures 
occurred.  The  progressive  failure  analysis  predicted  failures  around  the  hole  region  at 
the  stiffener  discontinuity.  The  final  failure  of  the  experiment  showed  that  local 
delaminations  and  debonds  were  present  near  the  hole,  and  edge  delaminations  were 
present  near  the  panel  midlength.  The  progressive  failure  analysis  results  did  not 
compare  well  to  the  test  results  since  delamination  failure  modes  were  not  included  in 
this  progressive  failure  approach. 

1  -2.2  Buckling  and  Postbuckling  Analvsis  of  Laminated  Composite  Plate  with 
Delamination 

1-2.2.1  One  dimensional  single  delamination 

Chai  et  al.  [43]  developed  a  one-dimensional  analytical  model  to  predict  through- 
the-width  delamination  buckling  and  growth  based  on  Euler  beam  theory.  Kardomateas 
[44]  investigated  effects  of  transverse  shear  and  end  fixity  of  delaminated  composite  by 
improving  the  model  used  Chai  et  al.  A  first-order  shear  deformation  theory  using 
variational  principle  was  proposed  by  Chen  [45].  Yin  and  his  group  [46-48]  investigated 
the  effects  of  bending-extension  coupling  on  postbuckling  behavior.  They  evaluated 
strain  energy  release  rate  by  using  J-integral  over  a  surface  that  encloses  the 
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delamination  boundary. 

Simitses  et  al.  [47]  developed  a  one-dimensional  model  similar  to  one  used  by 
Chai  et  al.  [42]  to  predict  critical  loads  for  delaminated  homogeneous  plates  with  both 
simply  supported  and  clamped  ends.  They  showed  that  the  buckling  loads  could  serve  in 
certain  cases  as  a  measure  of  the  load  carrying  capacity  of  the  delaminated 
configurations.  In  other  cases,  the  buckling  load  is  very  small  and  delamination  growth 
is  a  strong  possibility,  depending  on  the  toughness  of  the  material. 

Yin  et  al.  [48]  found  that  a  delamination  length  is  short  and  located  near  mid- 
plane  of  the  plate,  the  buckling  load  of  the  delaminated  plate  is  close  to  the  lower  bound 
of  the  ultimate  axial  load  capacity.  When  a  delamination  length  is  long  and  locates  near 
surface  of  plate,  the  postbuckling  axial  loads  can  be  considerably  greater  than  the 
buckling  loads,  while  the  failure  of  plate  may  or  may  not  be  governed  by  delamination 
growth. 

The  effects  of  bending-extension  coupling  as  well  as  imperfection  were 
investigated  by  Sheinman  and  Shouffer  [49].  They  found  that  the  coupling  effect 
reduces  the  load  carrying  capacity,  and  imperfection  sensitivity  of  global  postbuckling 
deformation  is  very  high. 

Wang  [50]  proposed  the  concept  of  a  continuous  analysis  for  determining 
interface  stresses  and  strain  energy  release  rate  for  the  delamination  at  the  interface  of 
skin  and  flange.  A  shear  deformable  beam  finite  element  with  nodes  offset  to  either  the 
top  or  bottom  side  was  proposed  by  Sankar  [51].  Kyoung  and  Kim  [52]  investigated 
asymmetric  delamination  with  respect  to  the  center  of  the  beam-plate.  In  their  study,  a 
variational  principle  based  on  shear  deformation  theory  was  used  to  calculate  buckling 
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load  of  orthotropic  laminated  beam-plate  with  through-the-width  delamination. 
Recently,  Gu  and  Chattopadhyay  [53]  carried  out  compression  tests  on  graphite/epoxy 
composites  plates  with  delaminations  to  evaluate  the  critical  load  and  the  actual 
postbuckling  load-carrying  capacity.  They  observed  that  composite  laminates  can  carry 
higher  loads  after  buckling.  For  the  particular  case  they  studied,  the  ultimate  load  is 
found  to  be  as  high  as  three  times  the  buckling  load. 
1.2.2.2  One-dimensional  multiple  delaminations 

Kutlu  and  Chang  [54-55]  investigated  the  compressive  response  of  composite 
laminates  with  multiple  delaminations.  They  found  that  multiple  delaminations  can 
reduce  the  load-carrying  capacity  more  compared  to  a  single  delamination.  Suemasu 
[56-57]  developed  closed  form  solution  for  linear  bifurcation  buckling  load  based  on 
energy  method.  He  found  that  in  the  case  of  multiple  delaminations,  size  and  location 
significantly  affect  the  buckling  load. 

Lee  et  al.  [58]  proposed  a  layer-wise  approach  for  computing  the  buckling  loads 
and  corresponding  buckling  mode  shapes.  It  was  found  that  the  anti-symmetric  buckling 
mode  is  dominant  for  a  composite  laminate  having  short  multiple  delaminations.  They 
also  addressed  the  effects  of  initial  imperfection  and  anisotropy  on  buckling  and 
postbuckling  response  of  delaminated  composite  plates.  An  analysis  procedure  for 
determining  the  buckling  load  of  beam-plates  having  multiple  delaminations  was  also 
presented  by  Wang  et  al.  [59]. 
1-2.2.3  Two-dimensional  single  delamination 

Whitcomb  [60]  studied  delamination  growth  caused  by  local  buckling  in 
composite  laminates  with  near  surface  delamination,  using  geometrically  nonlinear  finite 
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element  analysis.  He  showed  that  delamination  extension  does  not  occur  until  buckling 
is  significantly  progressed.  A  plane  finite  element  was  developed  by  Gim  [61]  based  on 
lamination  theory  that  included  the  effects  of  transverse  shear  deformation.  In  the 
modeling  of  two-dimensional  delaminations  in  laminated  plates,  the  undelaminated 
regions  was  modeled  by  a  single  layer  of  plate  elements  while  the  delaminated  region 
was  modeled  by  two  layers  of  plate  elements  with  node  offset. 

Sankar  and  Sonik  [62]  proposed  a  simple  expression  for  the  point-wise  strain 
energy  release  rate  along  the  delamination  front  using  Irwin's  crack  closure  technique. 
They  applied  this  technique  in  analyzing  stitched  double  cantilever  beam  specimens  as 
well  as  elliptic  delaminations  in  composite  plates. 

Klug  et  al.  [63]  investigated  efficient  modeling  of  postbuckling  delamination 
growth  using  plate  elements  and  gap  elements.  Energy  release  rate  was  computed  using 
virtual  crack  closure  technique.  From  this,  a  procedure  to  simulate  a  successive 
delamination  growth  was  proposed.  Kim  [64]  presented  a  modeling  approach  to  study 
the  postbuckling  behavior  of  composite  laminate  with  embedded  delamination  using 
two-dimensional  shell  element  and  rigid  beam  elements. 
1  -2.2.4  Two-dimensional  multiple  delaminations 

Suemasu  et  al.  [65-66]  analyzed  the  compressive  behavior  of  plates  with  mutiple 
delaminations  of  different  sizes.  They  showed  that  the  effect  of  variation  of  the  size  of 
delamination  on  the  compressive  behavior  is  significant  and  postbuckling  behavior  is 
different  from  that  of  plates  with  equal  sizes  of  delamination.  Zheng  and  Sun  [67] 
proposed  a  triple  plate  finite  element  model  to  analyze  delamination  interaction  in 
laminated  composite  structures.  Energy  release  rate  was  obtained  by  using  virtual  crack 
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closure  technique.  The  compatibility  conditions  between  interfaces  of  plates  were 
imposed  by  multi-point  constraint  equation.  The  results  for  delamination  interaction  in  a 
composite  laminated  circular  plate  under  three  point  bending  were  obtained. 

Lee  et  al.  [68]  developed  a  nonlinear  finite  element  code,  DELAM3D,  with  a 
three-dimensional  layered  solid  element  based  on  an  updated  Lagrangian  formulation. 
They  simulated  the  compressive  response  of  a  laminated  composite  plate  with  mutiple 
delaminations.  Contacts  of  delaminating  interfaces,  delamination  growth,  and  fiber- 
matrix  failure  were  also  considered  in  their  computation.  Double  cantilever  beam  (DCB) 
and  end  notched  flexure  (ENF)  tests  were  conducted  to  verify  the  energy  release  rate. 
Test  results  with  various  crack  numbers,  size,  location,  and  layer  orientation  compared 
well  with  the  numerical  results. 

1.3  Objective  and  Scope 

The  first  objective  of  this  study  is  to  develop  stiffened  panel  models  that  can  be 
used  to  predict  buckling  and  postbuckling  behavior  with  and  without  delamination.  The 
second  objective  is  to  investigate  the  delamination  growth  of  the  stiffened  panel  based 
on  fracture  mechanics  using  several  methods  of  computing  strain  energy  release  rate. 

Among  the  several  configurations  commonly  used  for  stiffened  panels  such  as 
hat-stiffened  panel,  J-stiffened  panel,  and  I-stiffened  panel,  blade  stiffened  panel  has  a 
simple  geometry  compared  with  other  stiffened  panels.  Therefore,  blade  stiffened  panel 
was  chosen  in  this  study.  However,  the  present  analysis  method  will  be  also  applicable 
to  any  kind  of  stiffener  configurations  that  have  a  flange  attached  to  the  skin. 
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Two  finite  element  modeling  approaches  for  stiffened  panels  are  commonly  used 
in  the  literature.  One  approach  is  to  model  the  skin  with  plate  elements  and  to  model  the 
stiffener  with  beam  elements.  The  other  is  to  model  both  skin  and  stiffener  with  plate 
elements.  The  second  approach  appears  more  attractive  for  modeling  the  debonded 
region  between  skin  and  flange.  In  this  study,  the  panel  skin  and  blade  stiffeners  are 
modeled  with  plate  elements.  The  nodal  penetration  of  the  delaminated  skin-stiffener 
interface  can  be  prevented  either  by  adjusting  spring  constants  of  fastener  elements  or  by 
gap  elements.  Furthermore,  an  energy  release  rate  for  calculating  delamination  extension 
is  computed  using  several  methods  based  on  fracture  mechanics. 

In  order  to  validate  the  present  modeling  approach,  a  plate  with  a  through-the- 
width  delamination  was  modeled  and  linear  bifurcation  and  nonlinear  postbuckling 
analyses  were  conducted.  Results  were  compared  with  the  available  experimental 
results. 

In  Chapter  2,  basic  finite  element  formulation  for  buckling  problem  and 
nonlinear  solution  algorithms  for  postbuckling  analysis  are  briefly  described.  Chapter  3 
provides  several  methods  for  computing  energy  release  rate  in  plate-like  structures  based 
on  fracture  mechanics.  In  Chapter  4,  effects  of  boundary  conditions,  material  properties, 
and  initial  geometric  imperfections  on  buckling  and  nonlinear  prebuckling  behavior  of 
blade  stiffened  panel  are  investigated.  Chapter  5  describes  modeling  of  delamination 
with  elastic  spring  element,  linear  buckling  analysis  results  of  both  delaminated 
composite  plates  and  debonded  stiffened  panels,  and  effects  of  delamination  length, 
stiffener  geometries  and  stacking  sequences  on  buckling  load.  Strain  energy  release  rate 
was  computed  using  the  strain  energy  derivative  method,  the  virtual  crack  closure 
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method,  and  the  crack-tip  force  method.  Finally,  Chapter  6  includes  the  summary  of  the 
present  study,  concluding  remarks,  and  suggestions  for  future  work. 


CHAPTER  2 

REVffiW  OF  THE  NONLINEAR  FINITE  ELEMENT  METHOD 
2.1  Nonlinear  Finite  Element  Formulations  Based  on  Continuum  Mechanics 

Stiffened  laminated  composite  panels  deform  continuously  under  compressive 
loads.  In  prebuckling  stage,  deformation  and  rotation  can  be  considered  as  infinitesimal 
in  general.  Thus,  prebuckling  response  of  stiffened  laminated  composite  panel  is  almost 
linear.  Classical  buckling  analysis  is  generally  used  to  estimate  the  critical  loads  of  stiff 
structures  such  as  the  Euler  column  subjected  to  axial  compression,  which  carry  design 
loads  by  axial  or  membrane  strength  rather  than  bending  strength.  The  out-of-the  plane 
deformation  before  buckling  is  therefore  almost  negligible  in  the  stiff  structures.  After 
buckling,  stiffened  laminated  composite  panels  exhibit  large  deformation  and  rotation. 
Thus,  nonlinear  formulation  is  required  in  order  to  include  the  effects  of  large 
deformation  and  rotation.  Many  researchers  [69-72]  have  efficiently  implemented 
general  nonlinear  finite  element  formulations  based  on  the  principles  of  continuum 
mechanics.  Two  different  approaches  have  been  used  in  incremental  non-linear  finite 
element  formulation.  The  first  approach  is  generally  called  Eulerian  or  updated 
formulation  where  static  and  kinematic  variables  are  referred  to  an  updated 
configuration  in  each  load  step.  The  second  approach  is  called  the  Lagrangian 
formulation,  where  all  static  and  kinematic  variables  are  referred  to  the  initial 
configuration.  The  updated  Lagrangian  is  more  suitable  for  analysis  of  the  structures 
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that  involve  very  large  deformations  while  the  total  Lagrangian  is  more  convenient  for 
analyzing  the  structures  with  moderately  large  deformations. 

The  primary  objective  of  this  section  is  to  review  briefly  the  non-linear  finite 
element  formulation  based  on  continuum  mechanics.  The  detailed  description  can  be 
found  in  References  [69,  70]. 

2.1.1  Equation  of  Equilibrium 

Using  the  principle  of  minimum  total  potential  energy  one  can  derive  the  finite 
element  equations.  Assume  that  there  exists  a  total  potential  energy  of  the  form  for 
linear  elastic  analysis 


where  {m}  is  the  displacement  vector,  {f^}  and  {/'}  are  body  and  surface  force  vectors. 
The  relationship  between  stresses  and  strain  is  of  the  form: 


where  [Q  is  the  constitutive  matrix.  The  condition  of  equilibrium  requires  that  the  first 
variation  of  the  total  potential  energy  vanish: 


n=-j{eV  {(j}dV -ijiuf  {f'}dV  +  j{uV  {f'}dS) 


(2.1) 


(2.2) 


^  =  \[S£nc7}dV-(j{Suf{f''}dV  +  l{Suf{f}dS)=0  (2.3) 


V  V 


From  Zienkiewicz  [69],  strain  can  be  expressed  in  matrix  notation  as 


{£}  =  [B]{u}  =  [B,]  +  [B,]{u} 


(2.4) 
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where  [Bo  ]  is  the  matrix  for  the  linear  infinitesimal  strain  and  matrix  [Bi  ]  contains  the 
nonlinear  strain  components. 

Using  Equation  (2.4)  we  can  rewrite  Equation  (2.3)  as: 

m  =  {Su}'{\[B]{a}dV  -  (/{/"      +  J{/'  }dS))  =  0  (2.5) 

V  V  s 

2. 1 .2  Eigenvalue  Buckling  Prediction 

The  stability  criterion  can  be  obtained  from  the  second  variation  of  the  total 
potential  energy.  If  the  second  variation  of  the  total  potential  energy  has  a  positive 
value,  then  a  system  is  stable.  Conversely,  If  the  second  variation  of  the  total  potential 
energy  has  a  negative  value  then  a  system  is  unstable.  Computing  the  second  variation 
of  total  potential  energy  from  Equation  (2.5)  as 

S'n  =  {Suf(lS[Bf{CT}dV  +  j[BYS{(jfdV)  (2.6.1) 
S'n  =  {Suf{lS[BY{a}dV  +  j[Bf[C][B]]d{ufdV)  (2.6.2) 

V 

From  Zienkiewicz  [69],  the  first  integral  of  Equation  (2.6.2)  can  generally  be  written  as 

lS[BY[(j}dV=[KJ{Su}  (2.7) 

where  [K„]  is  geometric  stiffness  matrix.  Substituting  Equation  (2.4)  into  the  second 

integral  of  Equation  (2.6.2)  and  rearranging,  the  second  variation  of  the  total  potential 
energy  can  be  written  as: 

S'n  =  {Suf[K^][Su}  (2.8) 
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In  Equation  (2.8),  the  tangential  stiffness  matrix  [Kt  ]  can  be  written  as 

[K,]^l[Bj[C][B,]dV  (2.9) 

V 

[K,]^j{[Bj[C][B,]  +  [Bj[C][B,]  +  [Bj[C][B,])dV  (2.10) 

V 

where  [Kq]  is  the  small  displacement  stiffness  matrix  and  [KJ  is  the  large 
displacement  stiffness  matrix. 

A  critical  point  is  obtained  when  the  tangent  stiffness  matrix  [Kj.]  has  at  least 
one  zero  eigenvalue.  The  stability  of  an  equilibrium  configuration  can  be  determined 
solving  the  eigenvalue  problem  at  the  current  equilibrium  state, 

[K^]{u'''}  =  ^}''{u''^}  (2.11) 
where  Z'''  is  the  rth  eigenvalue  and  {m^''}  is  the  corresponding  eigenvector. 

Computation  of  the  critical  point  must  be  done  in  two  steps.  First  the 
equilibrium  configuration  associated  with  a  given  load  level  P  is  computed.  Next  the 
stability  of  tangent  stiffness  matrix  is  examined  by  computing  the  eigenvalue  of  the 
tangent  stiffness  matrix  at  given  load  level  P.  This  method  of  determining  the  stability 
of  a  conservative  system  gives  accurate  results.  However,  it  is  computationally 
expensive  because  it  involves  the  solution  of  a  quadratic  eigenvalue  problem  for  the 
critical  load.  Linearized  buckling  analysis  calculates  critical  buckling  loads  based  on  a 
linear  extrapolation  of  the  structural  behavior  at  a  small  load  level.  Thus  it  is 
computationally  inexpensive.  From  the  fact  that  geometric  stiffness  matrix  [K^]  and 
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large  deformation  stiffness  matrix  [K^]  depend  on  the  load  level  P,  linearized  buckling 
analysis  approximates  the  tangent  stiffness  matrix  at  given  load  level  P  as  [40]: 

[K,{P)]  =  [K,]  +  ^i[K^{AP)]  +  [K,{AP)])  (2. 12) 

AP 

where  both  geometric  stiffness  matrix  and  large  displacement  matrix  are  computed  at 
small  load  level  AP .  If  we  assume  that  the  critical  load  is  equal  to  X-AP,  then  the 
condition  for  a  singular  point  becomes  a  standard  eigenvalue  problem.  There  are  two 
widely  used  numerical  methods  for  extracting  eigenvalues,  the  Lanczos  method  and  the 
subspace  iteration  method.  The  Lanczos  method  is  generally  faster  when  a  large 
number  of  eigenmodes  is  needed  for  a  system  with  many  degrees  of  freedom.  The 
subspace  iteration  method  is  effective  for  computing  a  small  number  of  eigenmodes. 

Based  on  the  assumption  that  the  displacements  {u}  are  infinitesimal  for  the 
small  load  AP ,  classical  buckling  problem  further  simplifies  Equation  (2.12)  as: 

([/^o]  +  ^'^'[^.(AP)]){m'^'}  =  0  (2.13) 
where  the  large  displacement  stiffness  matrix  [A:^]  in  Equation  (2.12)  is  ignored. 
However  the  applications  of  Equation  (2.13)  should  be  limited  in  practical  engineering 
problems.  In  order  to  avoid  the  erroneous  computation  of  the  stability  points  in  real 
engineering  applications,  the  stability  problem  should  be  investigated  using  full  tangent 
stiffness  matrix  in  Equation  (2.1 1)  [69]. 


2.2  Nonlinear  Solution  Methods 
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The  solution  of  nonlinear  finite  element  problems  includes  a  series  of  load  steps 
as  well  as  iterations  to  establish  equilibrium  at  the  new  load  level.  In  some  nonlinear 
static  analyses  the  equilibrium  configurations  corresponding  to  load  levels  can  be 
calculated  without  solving  for  other  equilibrium  configurations.  However,  when  the 
analysis  includes  path-dependent  nonlinear  conditions,  the  equilibrium  relation  needs  to 
be  solved  throughout  the  history  of  interest.  The  solution  may  be  obtained  by  using 
either  the  Newton-Raphson  or  the  modified  Newton-Raphson  methods.  The  Newton- 
Raphson  method  requires  evaluation  of  the  tangent  stiffness  matrix  at  each  iteration, 
which  is  computationally  expensive.  On  the  other  hand,  the  modified  Newton-Raphson 
method  evaluates  the  tangent  stiffness  matrix  at  each  load  step,  thus  improving  the 
computational  efficiency  compared  to  the  Newton-Raphson  method.  However,  the 
Newton  type  methods  fail  to  provide  a  solution  in  the  neighborhood  of  a  global 
bifurcation  or  limit  points  when  the  tangent  stiffness  matrix  becomes  singular  (see 
Figure  2.1). 

The  arc-length  method,  proposed  by  Riks  [73]  and  Wempner  [74],  and  modified 
by  Criesfield  [75,76],  is  an  effective  solution  procedure  to  search  equilibrium  path 
beyond  the  limit  points.  An  important  aspect  in  arc-length  methods  is  that  the  load  level 
is  treated  as  a  variable  in  addition  to  the  unknown  displacements  at  each  iteration  of  a 
load  step.  Thus,  an  additional  constraint  equation  comprising  the  displacements  and 
loads  is  required  to  calculate  the  load  level. 
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Figure  2. 1  Nonlinear  response  from  load  versus  displacement. 


2.2.1  Newton-Raphson  Method 


A  system  of  nonlinear  equilibrium  equation  can  be  written  as 

(m)  =/(«)-/ 
where  internal  forces  /(«)  is  defined  as 


I{u)  =  j[B]{(7}dV 


(2.14) 


(2.15) 


Unbalance  forces  4^(m)  represents  the  difference  between  internal  and  external  forces. 
The  basic  problem  is  to  find  solutions  that  satisfy  the  nonlinear  equilibrium  equation, 
^(u  )  =  0.  Since  Equation  (2.14)  cannot  be  solved  directly  for  the  displacement  of  u, 
both  an  incremental  equation  of  equilibrium  from  Equation  (2.14)  and  iterative 
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procedure  are  generally  used  for  its  solution.  The  Newton-Raphson  method  utilizes  the 
first-order  approximation  of  Equation  (2.14)  and  can  be  written  at  load  step  n+l  as 

4^  (<:',)  =  4^  =  0  (2.16) 

o  u 

Here  /  is  the  number  of  iteration,  and  the  tangent  stiffness  matrix  is  defined  as 

~^^^^^r  (2.17) 
ou  au 

From  Equation  (2.16)  we  have  the  following  iterative  correction  as 

(2.18) 

where  Kj  is  the  tangent  stiffness  matrix  at  the     iteration.  Thus  the  improved  solution 
can  be  computed  as 

"«Vi  =  u'n.^  +  Su[^,  (2.19) 
The  convergence  of  the  Newton-Raphson  method  is  generally  very  fast. 
However,  the  cost  of  computation  is  usually  high  due  to  the  calculation  and 
factorization  of  the  tangent  stiffness  matrix  at  each  iteration.  To  reduce  the  burden  of 
computational  cost,  the  modified  Newton-Raphson  was  introduced  in  which  the 
stiffness  matrix  is  approximated  as  a  constant:  K'j  ~  K,.   There  are  many  possible 

choices  of  the  approximated  stiffness  matrix  K, .  For  instance.  A",  can  be  chosen  as 
either  the  matrix  corresponding  to  the  first  iteration  or  some  previous  load  step. 
Schematics  of  Newton-Raphson  method  and  modified  Newton-Raphson  method  are 
shown  in  Figure  2.2. 
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Figure  2.2  Schematics  of  Newton-Raphson  and  modified  Newton-Raphson  methods. 
2.2.2  Arc -Length  Methods 

The  basic  feature  behind  the  standard  arc-length  method  is  that  load  level  A  is 
treated  as  a  variable  rather  than  constant  during  a  load  increment.  Thus  the  governing 
equilibrium  equation  (2.14)  can  be  rewritten  as 

=  /("„..)-  A„,,/o  =  0  (2.20) 

with 


(2.21) 


where  A  w  „  is  the  total  incremental  displacement  vector  at  the  n'^  iteration  and 
A  A,,  is  the  total  incremental  load  factor  at  iteration  n.  Since  the  load  level  is  treated 
as  a  variable,  we  need  an  extra  equation  that  constrains  the  iterative  displacements  to 
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follow  a  specified  path  towards  a  converged  solution.  The  constraint  equation 
proposed  by  Riks  [73]  may  be  expressed  as: 

(Am„)^(Am„)+  AU„/o7o  =  A/^  (2.22) 

where  A  /  is  a  user  defined  'incremental  length'  in  the  space  of  n+1  dimensions. 
Criesfield  [75]  has  proposed  a  modified  constraint  equation  that  includes  displacement 
components  alone  as: 

(Au„y  (Am  J  =  A/'  (2.23) 

It  is  possible  to  add  the  constraint  equation  (2.23)  to  the  system  of  equation 
directly  and  the  iterative  incremental  method  could  be  used  again,  however  this  may 
destroy  the  symmetry  and  banded  structure  of  the  equilibrium  equation.  Hence 
Criesfield  [75]  proposed  an  indirect  approaches  to  avoid  this  difficulty.  In  their 
approach,  the  displacements  at  a  given  iteration  /  is  written  as: 
Sul  =  -  K  ->       (A  A'„  +  SA 

=  -  ,-■  (T  (A  A'„  )  -  SA  (2.24) 
=  Su-:(AA'„)  +  SA  '„    (Su:),  dui  =  K  r  fo 

where  SiF^  are  the  iterative  displacements  corresponding  to  the  residual  forces  4^^^, , 
and  the  tangent  stiffness  matrix  K    '  is  formed  using  modified  Newton-Raphson  at 
the  beginning  of  each  increment  and  kept  fixed  for  all  iterations  within  the  increment. 
By  substituting  equation  (2.24)  in  which  SA  '„  is  still  undetermined  into  the  constraint 
equation  (2.23),  we  have 

(Am;-'  +  Su'JiAu:;'  +  Su')  =  Al'  (2.25) 


Expanding  Eq.  (2.25)  gives  a  quadratic  equation  for  the  unknown  iterative  load 
factor  SX    .  The  details  of  this  solution  procedure  are  given  in  Ref.  [75]. 
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.  Figure  2.3  One-dimensional  interpretation  of  spherical  arc-length  procedure. 

2.3  Solution  Strategy 

The  buckling  analysis  provides  information  about  the  load  level  at  which 
bifurcation  occurs.  In  some  cases  the  structure  withstand  far  above  the  buckling  load 
without  significant  damage.  In  other  cases  the  structure  collapses  well  below  this  load 
due  to  imperfection  sensitivity.  Stability  loss  at  a  bifurcation  point  occurs  only  if  the 
corresponding  deformation  mode  is  not  contained  in  the  deformation  mode  for 
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arbitrarily  small  loads.  For  example,  a  flat  plate  with  in-plane  loading  exhibits  no  lateral 
displacements  in  the  pre-buckling  range.  Thus  it  does  not  contain  the  lateral 
displacement  modes  of  the  buckled  plate.  Likewise,  if  the  structure  as  well  as  loading  is 
symmetric  about  some  plane,  all  deformation  modes  antisymmetric  with  respect  to  that 
plane  are  possible  bifurcation  buckling  modes.  In  such  cases,  we  may  choose  to 
perform  a  buckling  analysis  with  a  nonlinear  basic  stress  state.  Sometimes  when  a 
bifurcation  point  does  not  exist  at  all,  the  bifurcation  buckling  approach  may  still 
considered  as  an  acceptable  measure  of  the  critical  load.  If  the  structure  is  statically 
indeterminate  and  thus  allows  favorable  redistribution  of  the  stresses  (e.g.,  shells  with 
cutouts),  then  the  bifurcation  approach  is  too  conservative.  If  the  stiffness  of  the  shell 
deteriorates  with  increasing  load  (e.g.,  long  cylinders  under  bending),  the  bifurcation 
approach  gives  unconservative  results.  The  bifurcation  buckling  analysis  with  a  linear 
stress  state  is  probably  a  good  approximation  for  any  case  in  which  the  squares  of  the 
rotations  in  the  linear  solution  are  small  in  comparison  to  the  membrane  strains  at  the 
load  level  corresponding  to  bifurcation. 

Most  of  commercially  available  finite  element  codes  provide  an  option  to 
perform  the  stress  analysis  first,  save  the  data  for  a  certain  number  of  load  steps  on  tape 
and  later  decide  for  which  of  those  load  steps  buckling  loads  should  be  obtained  [IS- 
IS] This  option  may  save  some  computing  time.  First,  it  may  be  easier  to  decide  on  the 
load  levels  at  which  eigenvalues  are  desired  after  the  results  of  the  stress  analysis  have 
been  conducted.  Next,  it  makes  possible  to  find  additional  eigenvalues  in  a  subsequent 
run.  When  eigenvalues  are  computed  in  a  later  run,  the  data  deck  for  the  nonlinear 
prestress  analysis  can  be  used.  In  addition,  the  user  has  the  option  to  select  certain  data 


31 

sets  saved  on  tape  for  eigenvalue  solution.  Changes  in  boundary  conditions  are  also 
permitted  at  this  time.  In  most  practical  applications  one  range  of  eigenvalues  is 
particularly  important,  especially  to  a  sequence  of  the  lowest  eigenvalues. 

If  there  exists  a  symmetry  plane,  in  loading  as  well  as  in  geometry,  the  size  of 
the  problem  can  be  reduced  and  significant  savings  in  the  total  computational  effort  can 
be  achieved.  If  the  structure  on  one  side  of  the  symmetry  plane  is  considered,  only  the 
frequencies  of  symmetric  modes  are  obtained.  If  the  eigenvalue  analysis  is  based  on  a 
nonzero  prestress  analysis  with  symmetry  conditions  and  an  eigenvalue  analysis  with 
boundary  conditions  corresponding  to  anti-symmetry. 

The  eigenvalue  approach  for  bifurcation  buckling  analysis  with  linear  stress 
state  is  slightly  more  complicated  than  the  vibration  problem  because  eigenvalues  can 
be  negative  as  well  as  positive.  Often  the  analyst  is  only  interested  in  one  eigenvalue, 
the  lowest  positive  one.  If  the  analysis  is  performed  without  a  shift,  it  may  happen  that 
only  negative  eigenvalues  are  found  because  these  are  smaller  in  magnitude.  In  that 
case,  the  analysis  has  to  be  repeated  with  a  positive  shift.  In  choosing  the  shift  for  a 
repeated  run  the  user  can  utilize  the  fact  that  the  smallest  positive  eigenvalue  is  larger  in 
magnitude  than  the  largest  of  the  negative  eigenvalues  that  were  found.  Sometimes  the 
buckling  loads  are  symmetric  with  respect  to  zero.  This  is  the  case,  for  example,  if  a 
plate  or  a  cylinder  is  subjected  to  uniform  a  shear  load.  It  may  often  be  advisable  to 
request  more  than  one  eigenvalue  also  in  bucking  analysis.  If  the  structure  shows 
insufficient  strength  and  only  the  lowest  eigenvalue  and  corresponding  mode  are 
known,  reinforcements  may  be  introduced  that  have  little  effect  on  secondary  buckling 
mode  with  the  eigenvalue  below  the  design  load. 
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Nonlinear  analysis  is  computationally  expensive  compared  to  linear  analysis.  In 
order  to  get  a  sound  analysis  results  from  nonlinear  analysis,  analyst  should  have  better 
insight  into  the  behavior  of  the  analysis  model.  It  is  usually  possible  to  save  time  as 
well  as  computer  cost  by  preliminary  determination  of  approximate  values  for  the 
buckling  load,  under  negative  as  well  as  positive  load,  before  a  large  scale  analysis  is 
carried  out.  A  linear  analysis  with  a  rather  coarse  grid  will  give  some  idea  about  the 
stress  distribution  and  verify  the  nature  of  the  behavior  prior  to  executing  nonlinear 
analysis  with  refined  model.  The  size  of  the  finite  element  model  should  be  determined 
based  on  the  requirement  of  accuracy,  the  efficiency,  and  the  time  constraint.  Prior 
knowledge  of  the  geometric  modeling  will  increase  the  efficiency  of  an  analysis.  Type 
of  element  and  size  of  element  should  be  carefully  selected  to  obtain  high  accuracy. 
Further,  analyst  should  identify  the  type  of  nonlinearity  and  localize  the  nonlinear 
region  for  computational  efficiency.  To  identify  the  type  of  nonlinearity,  it  is  also 
helpful  to  examine  the  deformed  shapes  at  various  stages  of  loading  (pre-buckling, 
critical  buckling,  limit  load,  and  postbuckling)  [25]. 


CHAPTER  3 
COMPUTATION  OF  ENERGY  RELEASE  RATE 

3.1  Introduction 

Fracture  mechanics  concepts  have  been  successfully  applied  to  predict  the  loads 
which  initiates  the  delamination  extension,  and  also  for  predicting  their  stability.  The 
energy  release  rate  G  has  been  accepted  as  a  measure  for  predicting  delamination 
propagation.  In  the  context  of  fracture  mechanics,  the  delamination  extension  is  assumed 
to  occur  when  the  computed  G  is  greater  than  the  experimentally  determined  critical 
energy  release  rate  Gc-  Most  of  the  available  methods  of  computing  G  use  2-D  or  3-D 
solid  elements.  However,  It  is  computationally  very  expensive  to  use  solid  elements  for 
modeling  the  entire  complicated  structure  of  an  aircraft  or  an  automobile.  For  example, 
consider  2-D  plane  strain  elements  for  computing  stress  intensity  factor  to  estimate  the 
strain  energy  release  rate  of  a  double  cantilever  beam.  A  fine  mesh  must  be  used  around 
crack-tip  in  order  to  capture  the  stress  gradient  ahead  of  crack-tip.  Figure  3.1  shows  the 
amount  of  mesh  density  required  to  calculate  the  stresses  to  compute  the  stress  intensity 
factor.  Thus  a  simplified  beam,  plate  or  shell  theory  are  frequently  used  for  structural 
analysis  of  complicated  structures.  Therefore,  it  may  me  a  good  idea  to  use  the  same 
theory  to  calculate  the  energy  release  rate  also.  There  are  three  methods  commonly  used 
for  computing  energy  release  rate  using  2-D  or  3-D  solid  elements.  These  are  Strain 

33 


34 

Energy  Derivative  Method  (SEDM),  J-intergral,  and  Virtual  Crack  Closure  technique 
(VCCT).  SEDM,  first  proposed  by  Dixon  and  Pook  [77],  evaluates  the  change  of 
potential  energy  as  a  crack  progresses.  Implementation  of  SEDM  in  a  finite  element 
analysis  is  straightforward  [78].  However,  it  gives  only  an  average  value  of  the  energy 
release  rate  along  the  delamination  front.  Further,  this  method  requires  two 
computations  of  potential  energy,  before  and  after  crack  propagation.  A  direct  evaluation 
of  energy  release  rate  requiring  only  a  single  computation  was  proposed  by  Rice  [79]. 
This  involves  the  calculation  of  an  integral  on  an  arbitrary  path  surrounding  the  crack 
tip.  This  integral,  known  as  the  J-integral,  is  path  independent.  The  Virtual  Crack 
Closure  Technique  (VCCT),  proposed  by  Irwin  [80],  is  a  method  for  computing  energy 
release  rate  for  self-similar  crack  extension.  This  method  assumes  that  the  strain  energy 
release  during  the  crack  extension  is  equal  to  the  work  required  to  close  the  opened 
crack  surfaces.  Many  investigators  [51-68]  have  proposed  VCCT  for  computing  energy 
release  rate  using  beam  and  plate  elements.  Based  on  plate  theory,  Sankar  and  Sonik 
[62]  proposed  the  Point-wise  Strain  Energy  Density  Method  (PSEDM).  PSEDM 
suggests  that  the  point-wise  strain  energy  release  rate  along  the  crack  front  is  the 
difference  between  strain  energy  densities  behind  and  ahead  of  crack  front.  In  this 
chapter  a  new  method  called  Crack  Tip  Force  Method  (CTFM)  based  on  plate  theory  is 
introduced.  The  application  of  the  various  methods  of  computing  G  for  laminated 
composite  structures  is  discussed. 
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3.2  Strain  Energy  Derivative  Method 


The  strain  energy  derivative  method  utilizes  the  change  in  total  strain  energy,  U, 
with  change  in  crack  length  from  a  to  a  +Aa  (see  Figure  3.2).  The  energy  release  rate 
can  be  obtained  in  a  straight  forward  manner  for  the  case  of  displacement  control  as 


G  =  - 


dU 
dA 


^A 


const  —deflection 


(3.1) 


where  AA  is  the  increase  in  crack  area  due  to  change  Aa  in  crack  length. 
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As  shown  in  Equation  (3.1),  this  method  needs  two  analyses  to  compute  energy  release 
rate  and  Aa  should  be  small  enough  to  obtain  accurate  results  In  the  case  of  load  control, 
the  expression  for  G  is  modified  as 


dA  ^A 


(3.2) 

const .— force 


3.3  Path  Independent  J-Integral 

Consider  a  homogeneous  body  of  linear  or  nonlinear  elastic  material  without 
singularity  shown  in  Fig.  3.3.  The  strain  energy  density  Uo  is  defined  as 

e 

U,=U{e)  =  ja,jd£.j  (3.3) 

0 

where  a is  the  stress  tensor  and  e^j  is  the  infinitesimal  strain  tensor.  The  J- 
integral  to  compute  the  energy  release  rate  G  is  defined  as  [79] 

J  ^  l(U        -  (T,jnjU.^)ds    i  =  1,2  7  =  1,2  (3.4) 
r 

where  u.  is  the  displacement,  n,  are  the  direction  cosines  of  the  outward  normal  along 
the  path  T  .  The  indices,  /  and  j  or  x  and  z  are  used  interchangeably  for  convenience. 

Further,  summation  is  performed  over  repeated  indices.  An  application  of  divergence 
theorem  gives: 
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Differentiating  the  strain  energy  density, 


JL — (J  — y_ 


dx        d£;-  dx  dx 


(3.6) 


=  —cr,  [-— (-^)  +  -— (  J-)] 

2         dx    dxj        dx  dx. 

3    .3m,,        3    ,  3m,. 
=  (7  ■   ( — '-)  =  (cr,..  — '-) 

dxj    dx         dxj  dx 

The  area  integral  in  Equation  (3.5)  vanishes.  Therefore  Equation  (3.4)  is  equal  to  zero 
for  any  closed  contour  F  .  In  order  to  understand  the  zero  volume  integral  described  in 
subsequent  section  easily,  consider  the  conservation  integral  around  a  region  with 

singularity  as  shown  in  Fig.  3.4.  The  J  integral  along  closed  paths  Fj  through 
F4  surrounding  crack-tip  vanishes  as 


J  =  J(f/„n,  -  a.jnjU.^)ds  =  0 


(3.7) 


But  (T  ,^  =  0  and  n  ^  =  0  along  path      and  F^  .  Thus  the  integral  along 
F]  clockwise  and  the  integral  along  F3  counterclockwise  sum  to  zero. 

J  =  j(U        -  (7.jnjU.  Jds  +  \(U        -  G  ..n^u.  ^)ds  =  0  (3.8) 

From  equation  (3.8)  we  can  show  that  J  integral  along  path  F,  and  J  integral  along  path 
F3  have  the  same  value: 
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J  =  \      0^,-  O  ..n  .u.  ^)ds  =  -  J(t/„n,  -  G  .^n.u.  ^)ds 

(3.9) 


where  normal  vector  m  ^  is  in  the  opposite  direction  with  respect  to  normal  vector 


n  ^  in  path  F. 


n 


Figure  3.3  Conservation  integral  around  a  region  with  no  singularities. 


Figure  3.4  Conservation  integral  around  a  region  with  singularities. 
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3.4  Zero- Volume  J-Integral 


Consider  a  portion  of  the  delaminated  beam  as  shown  in  Fig.  3.5.  The  beam 
example  is  used  to  minimize  the  complexity  of  derivation  but  this  method  can  be 
extended  easily  to  delaminated  plates.  It  can  be  assumed  that  the  delamination  length  or 
crack  length  a  is  much  larger  that  the  thickness  h  of  the  thicker  sublaminate.  If  the  path 
of  the  integral  ABCDEF  is  away  from  crack-tip,  then  the  beam  theory  stresses  along  this 
path  are  reasonably  accurate  compared  to  exact  elasticity  solutions.  Further,  the  J- 
integral  will  vanish  along  the  two  horizontal  paths  BC  and  DE.  Thus  the  integral  is  given 
as  the  sum  of  integrals  along  the  three  vertical  paths:  AB,  CD  and  EF.  Next,  it  will  be 
shown  that  these  vertical  paths  can  be  moved  near  the  crack-tip  without  losing  any 
computational  accuracy  of  G. 

The  vanishing  of  the  J-integral  around  a  closed  path  in  an  elastic  material  under 
small  strain  assumptions  is  a  consequence  of  the  two  differential  equations  of 
equilibrium  satisfied  by  the  stress  components: 


dx 


XX 


+ 


dz 


dT 


zx 


=  0 


(3.10) 


+ 


zz_ 


=  0 


The  stress  field  in  a  laminated  beam  given  by  the  shear  deformable  beam  theory 
may  not  be  accurate  near  the  crack  tip,  however  the  stress  components  satisfy  the 
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equilibrium  equations  exactly.  This  is  because  that  the  transverse  shear  stresses  T  in 
the  beam  are  computed  actually  by  substituting  for  (7  ^  and  then  integrating  the  first 
equilibrium  equations.  Thus,  the  first  of  Equation  (3.10)  is  satisfied.  The  shear  stresses 
at  a  cross  section  are  proportional  to  the  shear  stress  T  ^,  ,  that  is  constant  along  each 
ligament  of  the  delaminated  beam  as  well  as  the  intact  beam  ahead  of  the  crack  tip.  Thus 
the  shear  stresses  T  are  independent  of  x  in  each  of  the  sublaminates  and  the  first 
term  in  the  second  equilibrium  equation  is  zero.  Since  beam  theory  assumes  that 
(T  ^  are  negligible,  the  second  term  is  also  zero.  Thus,  second  equilibrium  equation  is 
also  identically  satisfied.  Then  the  J-integral  evaluated  using  beam  theory  around  the 
closed  path  ABCDEF  in  Figure  3.5  is  identically  equal  to  zero.  Further  if  we  decompose 
delaminated  beam  with  three  sublaminates  as  shown  in  Figure  3.6  and  consider  J- 
integral  for  Sublaminate  2.  Because  the  integrals  along  the  horizontal  paths  are  zero,  it 

can  be  shown  that  Jab-Jhg-  Similarly  it  can  be  shown  that  J  cd  =  J  kl  and 
^  EF  =  J  MN  ■  Thus  it  is  now  possible  to  move  the  three  vertical  paths  AB,  CD  and  EF 
to  near  the  crack  tip  (HG,  KL  and  MN)  without  loss  of  accuracy.  The  J-integral 
evaluated  around  the  Paths  2,3,  and  1  (HGKLMN)  near  the  crack  tip  in  Figure  3.5  has 
been  called  the  zero-volume  J-integral  or  zero-area  J-integral,  and  it  is  given  by 

G  =  y"'  +  /^^>-Fy^^>  (3.11) 

where  superscript  (1),  (2),  and  (3)  represents  the  paths  I,  2,  and  3,  respectively. 
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Figure  3.6  Three  sublaminates  of  delaminated  beam. 
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3.5  G  from  Energy  Densities 

Consider  the  J-Integral  along  Path  1  in  Figure  3.5.  Along  this  path  n^^  -1  and 
n^=0.  Hence  the  Integral  can  be  written  as: 

N 


=  \{-U,+o^u^+T^w  Jds  (3.12) 


M 


We  will  add  and  subtract  to  the  integrand  in  Equation  (3.12): 

N  N 

J'''  =  \(-Uo+(^.u  ^+T^^{w        J)ds  -  jr^ujs  (3.13) 

M  M 

The  term  u,^  can  be  identified  as  the  rotation  at  the  crack  tip,  y/ ,  which  is  common  to  all 
the  three  paths  1,  2  and  3.  Further  the  second  and  third  terms  in  the  first  integral  in 
Equation  (3.13)  equal  to  2(7  g  Hence  Equation  (3. 13)  can  be  written  as: 

N  N 

y*"  =  \u,ds  -ii/,\T^ds  (3.14) 

M  M 

Equation  (3.14)  can  be  further  simplified  as: 

/''^  =  f//"-t^,y,  (3.15) 

where  U  and  V,  are  the  strain  energy  per  unit  length  and  shear  force  resultant  at  the 
cross  section  1.  Similar  results  can  be  derived  for  paths  2  and  3  as  follows: 


J">  =  -U,'"  +  y^,V, 

It  should  be  noted  that  in  Equation  (3.16)  switches  signs  because  of  change  in  the  sign  of 
Hx  from  -I  to  +1  for  the  Path  3.  Adding  all  the  three  integrals  and  nothing  that  the  shear 
force  resultants  must  satisfy  the  equilibrium  condition  V/+ V2  =  V'i ,  we  find  that 

=  f/l"+i/<^>-f/f'  (3.17) 

Thus  the  energy  release  rate  is  the  difference  between  the  strain  energy  densities  just 
behind  and  just  ahead  of  the  crack  tip.  The  strain  energy  density  in  the  context  of  beam 
refers  to  strain  energy  per  unit  length  of  the  beam,  U  ^  . 

3.6  G  in  terms  of  Crack  Tip  Force 

Consider  a  very  small  segment  of  the  beam  of  length  2  A  a:  surrounding  the 
crack-tip  (see  Figure  3.3).  It  will  be  convenient  to  shift  the  xz  coordinates  such  that  the 
xy  plane  coincides  with  the  plane  of  delamination.  Further,  we  will  divide  the  laminate 
into  4  sub-laminates  2  behind  and  2  ahead  of  the  crack  tip  as  shown  in  Figure  3.7.  Let 
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(3.16) 
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the  force  and  moment  resultants  near  the  crack  tip  in  any  sub-laminate  be  represented  by 
a  column  matrix  £  such  that  F_    =  [  P ,  M  ,  V  ]  where  P,  M  and  V  are  the  axial  force, 


Figure  3.7  Sub-laminates  in  a  delaminated  beam  and  the  coordinate  system. 

bending  moment  and  shear  force  resultants.  An  underscore  denotes  matrix  and  a 
superscript  T  denotes  transpose  of  a  matrix.  It  should  be  noted  that  the  force  and  moment 
resultants  are  resolved  about  the  x-axis,  which  lies  in  the  delamination  plane.  Thus  there 
is  an  offset  between  the  laminate  mid-planes  and  the  xy  plane.  The  force  resultants  in 
each  sublaminate  are  denoted  by  F/,  F2,  F5  and  F4.  The  compliance  matrix  of  the  top 

and  bottom  sub-laminates  will  be  denoted  by  C,  and  Q  .  The  deformation  in  a 
sublaminate  is  then  given  by 

£=C    F  (3.18) 

where  the  deformations  are  defined  by: 
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/  =K«''^.'r.J  (3.19) 
where  the  components  of  the  deformations  are  the  strain  along  the  x-axis  (not  the  sub- 
laminate  mid-plane),  rate  of  change  of  rotation  and  the  transverse  shear  strain, 
respectively.  The  force  resultants  are  related  by  the  equilibrium  conditions 

Ei_=  l±+  E±  (3.20) 
Further,  since  the  sub-laminates  3  and  4  are  intact  (not  delaminated)  the  deformations  in 
them  should  be  identical,  i.e.  ^3=^4,  and  hence 

£±El^£±E±  (3-21) 
If  F]  and  £2  are  given,  then  £5  and  F4  can  be  calculated  using  Equations.  (3.20),  and 
(3.21).  The  strain  energy  per  unit  length  in  any  sub-laminate  is  given  by: 


(3.22) 


L  2  

Substituting  Equation  (3.22)  into  Equation  (3.17)  we  obtain 

G  =  I  f  r  C ,  £r  +  i  F  [  C .  f ,  -  i  F I C .  £,  -  i  £:  C ,  F ,  (3.23) 

Using  the  relations  in  Equations  (3.20)  and  (3.21)  an  interesting  expression  for  G  can  be 
derived  as 

G  =  \{F_\-  F_\){C,  +  C,){F_,  -  F, )  (3.24) 

The  term  (£l4  -  £1, )  is  actually  the  force  transmitted  through  the  crack  tip  between  the 
top  and  bottom  sub-laminates,  and  can  be  called  the  crack-tip  force,  F , .  If  a  rigid  link 
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is  used  to  connect  the  top  and  bottom  crack  tip  nodes  in  a  finite  element  model,  then  the 
forces  transmitted  by  the  rigid  link  will  be  exactly  equal  to  the  above  crack-tip  forces.  It 
may  be  noted  that  the  crack  tip  force  vector  F^^  have  three  components,  an  axial  force,  a 
couple  and  a  transverse  force,  corresponding  to  each  degree  of  freedom  of  the  crack  tip 
nodes,  u,  y/  and  w. 

Another  important  implication  of  Equation  (3.23)  is  that  although  there  are  6 
independent  forces  P, ,  V, ,  M  , ,  P2  >  ^2  M  ^  that  can  be  applied  to  the  two 
delaminated  beam  ligaments  (see  Figure  3.5),  G  depends  only  on  three  crack  tip  force 
components  (see  Figure  3.8).  If  the  forces  F^i  and  2  such  that  =  e_2 ,  i.e., 
C_,E.\  =  C^h^i '  using  Equations  (3.18)  and  (3.19)  one  can  show  that  F_^  -  F_^,  and 
then  G  =  0  .  If  the  forces  on  the  top  and  bottom  sub-laminates  1  and  2  are  such  that 
they  produce  conforming  deformations  (£,  =£2)'  then  the  same  forces  act  in 
sublaminates  4  and  3,  respectively,  producing  conforming  deformations  ( £3  =  £4 ). 
Thus  there  is  no  need  for  any  interaction  between  the  top  and  bottom  laminates  at  the 
crack-tip,  and  hence  G  =  0  . 
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2-D  SOLID  BEAM 

MODEL 

Figure  3.8  Crack-tip  forces  in  beam 


3.7  Virtual  Crack  Closure  Technique 

The  virtual  crack  closure  technique  has  been  used  for  plate  and  beam  fracture  problems 
by  many  researchers.  We  will  derive  the  VCCT  from  the  Crack-Tip  Force  Method.  The 
expression  for  G  in  Equation  (3.23)  can  be  written  as: 

G=^Ll  (C ,  (F ,  -  F , )  -f-  C ,  (F  3  -  F  3 ))  (3.25) 

where  F,c  is  the  matrix  of  crack  tip  forces,  and  Equation  (3.18)  is  used  in  deriving 
Equation  (3.25).  Using  the  compatibility  quation  (3.19)  in  (3.24),  we  obtain 

G  =  ^FJA-C,F_,  +C,F:,)  (3.26) 


Since  C£.denote  the  deformations  we  can  write  (3.26)  as 


49 


2' 


(2) 


,(2) 


(3.27) 


In  deriving  the  last  term  of  the  column  matrix  in  Equation  (3.27),  we  have  used  the  fact 
that  the  beam  rotation  at  the  crack  tip  is  same  for  both  ligaments  1  and  2,  i.e.  ^,  =  y/j- 
Multiplying  and  dividing  the  right  hand  side  of  Equation  (3.27)  by  -  Ax  where  Ax  is 
a  small  length  used  in  the  virtual  crack  closure  method,  we  obtain 


G  = 


lAx' 


(«;'^-«j'^)-(uj^>-«5'>) 

(^(»_^('))_(^(2)_^(0) 
(W0)+^(0)_(^(2)_^,0) 


(3.28) 


The  superscript  {t)  in  Equation  (3.28)  denotes  displacements  and  rotation  at  the  crack 
tip,  and  superscript  (/)  and  (2)  denote  respectively  the  displacements  of  the  top  and 
bottom  ligaments  at  a  distance  Ax  from  the  crack  tip.  In  deriving  Equation  (3.28)  we 
have  used  the  finite  difference  approximation  of  the  type. 
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Canceling  the  crack  tip  displacements  in  Equation  (3.28)  we  obtain  the  equations  for  the 
virtual  crack  closure  method  as: 


G  = 


2Ax 


(3.30) 


3.8  Extension  to  Delaminated  Plates 

In  the  case  of  delaminations  in  a  plate  the  energy  release  rate  G  varies  along  the 
delamination  front.  Formulas  similar  to  Equations  (3.10)  (SEDM)  and  (3.22)  (VCCT) 
were  derived  by  Sankar  and  Sonik  [62].  In  this  section  we  will  derive  an  additional 
result  for  G{s)  similar  to  Equation  (3.23)  derived  for  beam.  We  can  use  the  same 
notation  as  we  used  for  beams  with  the  understanding  that  there  are  eight  force  and 
moment  resultants  and  eight  deformation  components: 


[F']  =  [N^  Q^] 
[^']  =  [^.o    ^,0    r^o  S  7^  ] 


The  laminate  compliance  matrix  [Q  will  be  an  8x8  symmetric  matrix,  and  it  relates  the 
force  resultants  and  deformations: 


51 

e  =  C_F_  (3.32) 


A 

B 

0 

B 

D 

0 

{F} 

0 

0 

K 

where,  the  [A],  [B]  and  [D]  are  the  classical  3x3  laminate  stiffness  matrices  and  [K]  is 
the  2  by  2  transverse  shear  stiffness  matrix.  In  the  context  of  plates  the  strain  energy 
density  is  defined  as  strain  energy  per  unit  area  of  the  plate  and  is  given  by 

U,=^fJcf:  (3.33) 

A  formula  for  G(s)  similar  to  that  in  Equation  (3.23)  is  given  by 

G(s)  =  (Ul!\s)  +  U\'\s)-U['\s)-U\'\s))  (3.34) 
where  the  superscripts  denote  the  four  sub-laminates  behind  and  ahead  of  the 
delamination  front,  and  s  denotes  the  location  of  the  point  along  the  delamination  front. 
Sub-laminates  1  and  4  are  above  the  delamination  plane,  and  2  and  3  are  below  the 
delamination  plane.  From  Equation  (3.34)  one  can  derive  another  expression  for  G(s)  as: 

G  =  y  (£.:  -  Fj  ){C,  +  C,  )(F ,  -  F ,  )  (3.35) 
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The  derivation  of  Equation  (3.35)  is  very  similar  to  that  of  Equation  (3.23).  As  before, 

the  term  ( -  E-i)  is  the  matrix  of  crack-tip  forces.  They  also  represent  the  jump  in 
force  and  moment  resultants  that  occur  across  the  delamination  front. 

Sankar  and  Sonik  [62]  showed  that  three  of  the  eight  force  resultants  will  be 
continuous  across  the  delamination  front.  Assume  a  coordinate  system  such  that  the  x- 
axis  is  normal  to  the  crack  front,  y-axis  is  tangential  to  the  crack  front  and  z-  is  the 

thickness  direction.  Then  the  continuous  force  resultants  are  N  ^  and  Q  ^ .  Thus 
the  jumps  in  these  force  resultants  are  zero,  i.e.. 


^(4)_^(.)^^(2)_^(3)^Q 

y  y  y  y 

Ml" -Ml"  =M'-;> -Ml"  =0 


Thus  there  will  be  only  five  components  to  the  crack  tip  forces  (see  Figure  3.9):  three 
forces  in  the  x,  y  and  z  directions;  two  couples,  about  x  and  y  axes,  respectively.  The 
three  forces  will  be  the  jump  in  !^  ^  and  across  the  delamination  front,  either 
in  the  top  laminates  (1  and  4)  or  bottom  laminates  (2  and  3).  The  two  crack  tip  couples 
are  the  jumps  in  M  ^.  Since  the  jumps  in  N  ^  and  Q  ^  are  equal  to  zero  and 
they  do  not  contribute  to  the  crack-tip  forces,  we  can  delete  the  2"<i,  S^h  and  7^^  j-o^s 
and  columns  in  Q  and  C^;  we  will  denote  them  by  C]  and  C  'y. 
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J 


Figure  3.9  Crack-tip  forces  in  delaminated  plate. 


Then  from  Equation  (3.33)  an  expression  of  point-wise  energy  release  rate  can  be 
derived  as: 


where  the  crack  tip  forces  F^.  are  given  by: 

[Fl  ]  =  [( ^  -       )^  (N^:^  -  < ),  (M r ^  -  M ), 


(3.37) 


(3.38) 


The  compliance  matrices  C  will  take  the  form: 


c 

Ml 

M3 

M4 

M6 

Ms 

c 

'-'13 

c 

c 

M6 

Ms 

c 

Q4 

Q6 

Ms 

c 

M6 

Q 

Q6 

c 

'-68 

c 

c 

Ms 

'-'48 

'-68 

Ms 

(3.39) 


where  the  Cij  are  the  coefficents  of  the  full  compliance  matrix  C,  or  .  It  should  be 
mentioned  that  the  xyz  coordinates  should  be  moved  along  the  delamination  front  while 
using  Equation(3.30)  for  computing  pointwise  G(s). 

We  have  derived  three  formulars  Equations  (3.27),  (3.31),  and  (3.34)  for 
computing  the  point-wise  energy  release  rate  in  a  delaminated  plate.  Out  of  the  three. 
Equations  (3.32)  and  (3.35)  are  exact,  and  their  accuracy  is  limited  only  by  the  methods 
used  compute  the  force  and  moment  resultants  or  the  strain  energy  densities  ahead  and 
behind  the  delamination  front.  The  accuracy  of  the  VCCT  given  by  Equation  (3.30)  is 
limited  by  not  only  the  crack  tip  forces  but  also  the  mesh  size  which  will  define  the 
length  of  virtual  crack  growth. 


3.9  Summary 

In  this  chapter  a  new  method  called  Crack  Tip  Force  Method  (CTFM)  was 
introduced  and  derived.  Three  methods  of  computing  G  for  laminated  composites 
structures  are  also  discussed.  A  Crack  Tip  Force  Method  is  derived  for  computing 
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point-wise  energy  release  rate  along  the  delamination  front  in  delaminated  plates. 
Actually  the  method  can  be  derived  from  the  Virtual  Crack  Closure  Technique  or  the 
previously  derived  Strain  Energy  Density  Method.  However  the  CTFM  is 
computationally  simple  as  the  G  is  computed  using  the  forces  transmitted  at  the  crack- 
tip  between  the  top  and  bottom  sub-laminates  and  the  sub-laminate  properties.  An 
evaluation  of  the  aforementioned  methods,  their  applicability  to  general  laminates 
containing  delaminations,  and  debonded  stiffened  panels  will  be  presented  in  chapter  5. 


CHAPTER  4 

ANALYTICAL  AND  EXPERIMENTAL  CORRELATION  OF  A  STIFEENED 
COMPOSITE  PANEL  IN  AXL\L  COMPRESSION 

4.1  Introduction 

Buckling  and  imperfection  sensitivity  are  expensive  to  calculate  with  general 
finite  element  models.  Consequently,  the  optimization  of  stiffened  panels  often  employs 
simplified  models  that  are  exact  only  for  idealized  geometries  and  boundary  conditions 
(e.g.,  PASCO  [10],  or  PANDA2  [15-17]). 

Nagendra  et  al.  [33]  studied  the  optimum  design  of  blade  stiffened  panels  with 
holes  subjected  to  buckling  and  strain  constraints.  They  used  the  panel  analysis  and  sizing 
code  (PASCO),  based  on  a  linked  plate  model,  for  the  buckling  analysis  and  optimization 
with  continuous  thickness  design  variables,  and  the  Engineering  Analysis  Language 
(EAL  [34])  finite  element  analysis  code  for  calculating  strains  and  their  derivatives  with 
respect  to  design  variables.  Later,  the  optimally  designed  panels  with  and  without 
centrally  located  holes  were  tested,  and  analytical  and  experimental  results  were 
compared  [35].  Nagendra  et  al.  [80]  continued  the  optimum  design  study  of  blade 
stiffened  panel  using  PASCO  for  analysis,  and  a  genetic  algorithm  (GA)  for  the 
optimization  of  the  panel  laminate  stacking  sequences.  Several  designs  obtained  with  GA 
were  about  8%  lower  in  weight  compared  to  designs  previously  obtained  with  a 
continuous  optimization  procedure. 
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Recently,  three  of  the  panels  designed  by  Nagendra  et  al.  were  fabricated  and 
tested  by  the  Structural  Mechanics  Branch  at  NASA  Langley  Research  Center.  The 
experimental  failure  loads  differed  by  up  to  about  10%  from  the  design  load.  However, 
there  were  significant  differences  in  loading  and  boundary  conditions  between  the  design 
conditions  and  the  test  conditions. 

The  principal  objective  of  this  chapter  is  to  understand  the  effects  of  differences 
between  the  simplified  assumptions  made  in  the  design  model  and  the  actual  test 
conditions.  Another  objective  is  to  assess  the  effectiveness  of  the  simplified  PASCO 
model  originally  used  to  design  the  panel,  and  this  is  aided  by  comparing  the  results 
obtained  from  several  structural  analysis  programs  including  PANDA2,  STAGS,  and 
ABAQUS. 


4.2  Stiffened  Panel  Definition 

The  basic  configuration  of  the  panel  -  designated  as  the  baseline  design  - 
corresponds  to  the  design  in  the  9^^  row  in  Table  7  of  Reference  [80].  This  panel, 
designated  as  GA2461  (referring  to  the  design  weight  of  24.61  lb.),  is  30-inches  long  and 
32-inches  wide  with  four  equally  spaced  blade  stiffeners  (see  Figure  4.1).  The  laminates 
used  in  Reference  [95]  for  the  skin,  stiffener  blade,  and  stiffener  flange  for  the  baseline 
design  were  balanced,  symmetric  laminates  consisting  of  0°,  ±45°  and  90°  plies.  The 
skin  has  40  plies  with  a  stacking  sequence  [±45/90/1453/907+45,]^  and  the  stiffener 
flange  and  blade  have  an  identical  stacking  sequence  of  [±452/(±45/04)2 
/902/04/(±45/02)2/02/±45]s  with  a  total  of  68  plies.  Properties  of  the  AS4/3502  graphite 
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epoxy  material  used  in  Reference  5  are  given  in  Table  4. 1 . 

Table  4.1  Hercules  AS4/3502  graphite  epoxy  lamina  material  properties. 


Young's  modulus 
(longitudinal) 

E,  =18.50  X  10^  psi 

Young's  modulus 
(transverse) 

El  =1.64  X  lO^psi 

Shear  modulus 

G,2  =0.87  X  10^  psi 

Poisson's  ratio 

Vi2  =  0.3 

Density 

p  =0.057  lb  in-^ 

Ply  thickness 

tpiy  =  0.0052  in 

The  baseline  design  was  designed  to  support  an  axial  load  N  of  20,000  Ib./in.  In 
addition,  in  order  to  account  for  off  design  conditions,  imperfections  and  modeling 
inaccuracies,  a  shear  load  (N  ^  =  5000  Ib./in)  and  a  longitudinal  bow  type  (3%  of  the 
panel  length)  imperfections  were  added.  The  baseline  design  panel  was  assumed  to  be 
simply  supported  along  the  four  edges,  which  is  the  only  boundary  condition  that  can  be 
accurately  modeled  using  PASCO. 
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Figure  4. 1  Blade  stiffened  panel  with  four  equally  spaced  stiffeners  under  compression 
and  shear  loads.  All  dimensions  are  in  inches. 
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4.3  Test  Specimen  and  Test  Procedures 

The  test  specimens  were  fabricated  from  Hercules  AS4/3502  unidirectional 
graphite/epoxy  preimpregnated  tape  material.  The  skin  (32-in.  x  32-in.)  and  the  stiffeners 
were  cured  separately  in  an  autoclave.  The  stiffeners  were  machined  to  a  length  of  32 
inches,  and  then  bonded  to  the  skin  with  FM-300  film  adhesive.  The  panel  edges 
perpendicular  to  the  stiffeners  were  potted  with  an  aluminum  filled  epoxy  resin  to  prevent 
end  failure.  The  length  of  the  potted  area  was  1  in.  on  each  side.  Thus  the  effective  gage 
length  of  the  test  specimen  was  reduced  to  30  inches. 

The  test  specimen  was  loaded  in  compression  using  a  1 ,000,000-lb-capacity 
hydraulic  testing  machine.  The  specimen  was  flat-end  tested  without  lateral  edge 
supports,  and  no  deliberate  imperfection  was  introduced.  Electrical  resistance  strain  gages 
were  used  to  monitor  the  strains,  and  direct  current  differential  transformers  (DCDTs) 
were  used  to  monitor  longitudinal  in-plane  and  out-of-plane  displacements  at  selected 
locations  as  shown  in  Fig.  4.2.  All  electrical  signals  and  corresponding  applied  loads 
were  recorded  automatically  at  regular  time  intervals  during  the  tests. 

4.4  Linear  Buckling  Analysis 

4.4.1  PANDA2  and  STAGS 

The  analyses  performed  in  this  study  include  both  buckling  and  nonlinear 
postbuckling  calculations.  Linear  buckling  analyses  for  the  baseline  design  were 
conducted  by  using  both  PANDA2  and  STAGS.  Input  files  for  STAGS  linear  buckling 
analysis  were  generated  by  PANDA2.  Next,  the  effect  of  the  shear  load  and  the 
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DCDTs  (1-7)  for  out-of-plane  deformation  of  skin. 
DCDTs  (8-11)  for  out-of-plane  deformation  of  stiffener. 
DCDT  (12)  for  end-shortening  displacement 
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Figure  4.2  Layout  of  the  displacement  measurement  instrumentation  for  the  test  panel. 


Layer  without 
stiffness 


Reference  plane  (ABAQUS) 


0  0  0  0  Q-Q-e-o-- 


Reference  plane  (STAGS) 


Figure  4.3  Reference  plane  of  ABAQUS  and  STAGS  model. 
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imperfection  on  the  buckling  loads  of  the  baseline  design  were  investigated  using 
PANDA2,  which  employs  analysis  techniques  with  similar  level  of  fidelity  to  that  of 
PASCO.  In  PANDA2,  local  and  general  buckling  loads  are  calculated  by  either  closed- 
form  expressions  or  by  discretized  models  of  panel  cross  sections  based  on  an  energy 
method  [15]. 

STAGS  is  a  finite  element  code  for  general  purpose  analysis  of  shell  structures 
of  arbitrary  shape  and  complexity  [25].  STAGS  has  a  variety  of  finite  elements  suitable 
for  the  analysis  of  stiffened  plates  and  shells.  Four  node  quadrilateral  plate  elements 
with  cubic  lateral  displacement  variations  (called  410-  and  41 1-Elements)  are  efficient 
for  the  prediction  of  buckling  response  of  thin  shells.  For  thick  plates  in  which 
transverse  shear  deformation  is  important,  the  assumed  natural  strain  (ANS)  nine  node 
element  (480-Element)  can  be  selected  [16].  The  panel  investigated  here  warrants  the 
use  of  480-Element,  however  41 1 -Element  was  also  used  as  the  panel  was  designed  by 
PASCO,  which  does  not  model  shear  deformation.  STAGS  results  were  post-processed 
by  PATRAN,  which  is  a  commercial  software  for  pre-  and  post-processing  of  finite 
element  simulations  [82]. 

4.4.2  Finite  Element  Model 

The  STAGS  finite  element  model  for  the  panel  had  a  total  of  20  branched  shell 
units,  and  each  branched  shell  unit  had  65  x5  nodes  (for  the  32-in.  long  panel)  or  61  x  5 
nodes  (for  the  30-in.  long  baseline  design  panel),  respectively.  The  axial  compressive 
design  load  (640,000  lb)  was  applied  with  a  uniform  end-shortening  constraint  along 
with  compatibility  conditions  for  adjacent  shell  unit  interfaces.  In  the  test,  the  load  was 
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introduced  through  the  potted  ends.  To  simulate  this  boundary  condition,  the 
displacement  along  the  z-direction  and  the  rotation  along  the  x-direction  were 
constrained  at  the  nodal  points  in  the  potted  region.  The  adhesive  film  used  to  bond  the 
stiffener  to  the  skin  in  the  test  panel  was  modeled  by  adding  an  isotropic  layer  to  the 
model  to  simulate  the  bondline  between  the  skin  and  flange  with  a  thickness  of  0.0121 
inches.  The  skin  middle  surface  was  used  as  reference  surface  on  which  the  nodes  lie, 
and  the  offset  distance  of  the  middle  surfaces  of  the  skin-flange  combination  was 
modeled  as  an  eccentricity. 

An  alternative  finite  element  modeling  approach  with  ABAQUS  suggested  by 
Greene  of  HKS,  Inc.  was  also  used.  In  this  method  instead  of  locating  the  reference 
plane  at  the  mid-plane  of  skin,  the  bottom  plane  of  blade  stiffener  was  taken  as  the 
reference  plane.  In  order  to  handle  the  offset  distance  of  the  mid-plane  of  skin,  as  well 
as  skin-flange  combination,  an  additional  0°  ply,  with  negligibly  low  stiffness  was  added 
to  both  skin  and  skin-flange  laminates,  as  shown  in  Figure  4.3.  Both  the  nine-node  thin 
shell  element  (S9R5)  and  the  four-node  general  shell  element  (S4)  in  ABAQUS  were 
selected  for  the  stiffened  panel  models.  Both  shell  elements  can  account  for  transverse 
shear  deformations  and  should  therefore  be  compared  to  the  480-Element  in  STAGS. 
Instead  of  applying  a  compressive  load  at  the  panel  end,  uniform  compressive 
displacements  were  applied  at  the  nodal  points  along  the  loaded  edge.  In  order  to  ensure 
a  uniform  state  of  stress  along  the  entire  panel  length  and  also  to  prevent  bending  during 
the  pre-buckling  stage,  the  incremental  boundary  condition  option  available  in 
ABAQUS  was  chosen.  The  buckling  load  factor  was  computed  from  the  sum  of  the 
reaction  forces  at  the  boundary  node  set. 
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4.5  Results  of  Linear  Buckling  Analysis 

In  this  section  the  effects  of  geometric  imperfections,  additional  in-plane  shear 
loads,  boundary  conditions  and  material  property  variations  on  the  buckling  load  of  the 
stiffened  panel  are  discussed.  The  results  are  intended  to  shed  light  on  probable  reasons 
for  the  discrepancies  between  predicted  buckling  loads  and  corresponding  experimental 
results.  Furthermore,  the  effect  of  the  assumed  imperfections  and  the  addition  of  in- 
plane  shear  loads  on  the  robustness  of  the  design  is  also  noted. 

4.5.1  Effect  of  Geometric  Imperfections  and  Shear  Load 

A  summary  of  the  local  buckling  load  factors  with  and  without  shear  load,  and 
with  and  without  the  initial  bow  type  geometric  imperfections  (±3  %  of  the  panel  length) 
obtained  from  PANDA2  is  given  in  Table  4.2.  The  first  row  in  Table  4.2  also  includes  a 
comparison  of  PANDA2  results  and  the  STAGS  results  (both  480-  and  4 11 -Elements) 
for  the  perfect  panel  without  the  shear  load.  It  may  be  noted  that  the  PANDA2  results, 
both  Koiter  type  analysis  and  B0S0R4*  analysis,  agree  well  with  the  STAGS  411- 
Element  results.  The  buckling  mode  for  the  perfect  panel  obtained  using  STAGS  480- 
Element  is  shown  in  Figure  4.4.  According  to  PANDA2  results,  the  lowest  buckling  load 
corresponded  to  local  buckling,  which  suggests  that  the  differences  in  boundary 
conditions  between  the  analysis  and  the  experiment  will  not  have  a  large  effect  on  the 
results.  The  10%  difference  between  the  results  for  the  480-Element  and  the  411- 

*  B0S0R4  analysis  routine  in  PANDA2  calculates  local  buckling  load  for  the  single  panel  module  from  B0S0R4- 
type  strip  theory.' 
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Element  (Table  4.2)  is  probably  due  to  transverse  shear  deformation  since  the  thickness 
of  the  skin-flange  combination  is  0.56  inch.  Shear  deformation  was  not  included  in  the 
original  panel  design,  and  this  difference  indicates  that  that  effect  is  substantial. 

The  panel  with  negative  bow-type  imperfection  had  a  concave  surface  in  the 
middle  of  the  panel.  Thus,  the  blade  tip  is  subjected  to  less  axial  compression  and  skin  is 
under  more  axial  compression  than  that  of  the  perfect  panels  in  the  neighborhood  of 
mid-length  in  the  axial  direction.  Similarly,  the  blade  tip  near  the  boundary  is  under 
more  compression  and  the  skin  near  the  boundary  is  under  less  compression  than  that  of 
the  perfect  panel.  The  opposite  holds  for  the  positive  bow-type  imperfection.  From 
Table  4.2,  it  is  clear  that  the  effect  of  the  shear  load  on  the  buckling  load  is  small,  but  the 
effect  of  the  imperfection  is  very  significant.  From  the  last  two  rows  of  Table  4.2  one 
can  note  that  a  3%  positive  imperfection  results  in  a  very  low  buckling  load  factor.  The 
buckling  load  factor  reduces  from  1.256  to  0.394.  A  3%  negative  imperfection  also 
reduces  the  buckling  load  (from  1.256  to  0.856),  but  the  reduction  is  smaller  than  that 
for  a  positive  imperfection.  It  should  be  mentioned  that  a  3%  imperfection  is  very  large 
for  a  32-in.  stiffened  panel,  and  thus  will  lead  to  very  conservative  designs. 

4.5.2  Effects  of  Boundary  Condition  and  Material  Properties 

There  were  slight  differences  in  material  properties,  panel  dimensions  and  the 
boundary  conditions  between  the  baseline  deign  and  the  actual  test  conditions.  In  order 
to  understand  the  effects  of  these  differences,  analyses  were  carried  out  using  both  sets 
of  input  data.  The  differences  in  material  properties  and  dimensions  are  summarized  in 
Table  4.3  and  Table  4.4,  respectively.  While  the  changes  in  the  material  properties  can 
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be  input  directly,  the  differences  in  the  thickness  of  the  skin  or  flange  are  accounted  by 
implementing  a  proportional  change  in  the  ply  thickness  in  the  model.  A  detailed 
discussion  of  this  procedure  can  be  found  in  Ref.  4.  The  results  in  Table  4.5  indicate  that 
the  effects  of  boundary  conditions  and  material  properties  on  the  buckling  load  factors 
are  not  very  significant.  Comparison  of  the  first  two  rows  of  Table  4.5  reveals  the  effects 
of  changes  in  the  boundary  conditions.  Similarly,  results  in  the  last  two  rows  show  the 
effects  of  changes  in  material  properties  and  panel  dimensions.  The  buckling  mode 
shape  of  the  baseline  design  (simply  supported  on  4  sides)  predicted  by  STAGS  is 
shown  in  Figure  4.4.  The  overall  buckling  mode  shape  obtained  from  ABAQUS  (Figure 
4.5)  agrees  well  with  that  of  STAGS.  The  computed  lowest  buckling  load  factor  is 
slightly  higher  than  that  of  STAGS  (1.218  for  ABAQUS  vs.  1.168  for  STAGS).  This 
small  difference  may  be  due  to  modeling  differences  as  discussed  in  the  following 
section.  The  STAGS  prediction  of  the  buckling  mode  shape  of  the  test  panel  with  potted 
boundary  conditions  (other  two  edges  being  free)  is  shown  in  Figure  4.6. 


Table  4.2  Summary  of  the  local  buckling  load  factor  from  PANDA2  and  the  lowest 
buckling  load  factor  from  STAGS  (4  edges  simple  supported). 


PANDA2 
(Koiter  analysis) 

PANDA2 
(B0S0R4  analysis) 

STAGS 

(480 
element*) 

STAGS 

(411 
element) 

Loading  combination  with/without 
imperfection 

Panel  end 

Panel 
mid-length. 

Panel 
end 

Panel 
mid-length. 

N,=20,000  lb/in,  N,y=0  without 
imperfection 

1.256 

1.256 

1.328 

1.328 

1.168 

1.296 

N,=20,0001b/in,  N,y=5000  lb/in  without 
imperfection 

1.234 

1.234 

1.048 

1.048 

N.=20,0001b/in,  N,y=5000  lb/in  with 
imperfection  (+3%) 

1.234 

0.356 

1.048 

0.346 

N,=20,0001b/in,  N,y=5000  lb/in  with 
imperfection  (-3%) 

1.234 

0.856 

1.048 

0.920 

N,=20,000lb/in,  N,y=0  lb/in  with 
imperfection  (+3%) 

1.256 

0.394 

1.328 

0.398 

N,=20,0001b/in,  N.y=0  lb/in  with 
imperfection  (-3%) 

1.256 

1.026 

1.328 

1.083 

*includes  shear  deformation. 


Table  4.3  Material  properties  of  baseline  designs  and  test  specimen. 


E„(Msi) 

E,,(Msi) 

G„  (Msi) 

Vl2 

V21 

Skin 

Design 

18.50 

1.64 

0.87 

0.3 

0.0270 

Test 
specimen 

17.333 

1.64 

0.8151 

0.3 

0.0284 

Blade 

Design 

18.50 

1.64 

0.87 

0.3 

0.0270 

Test 
specimen 

19.125 

1.64 

0.8994 

0.3 

0.0257 

Flange 

Design 

18.50 

1.64 

0.87 

0.3 

0.0270 

Test 
specimen 

19.593 

1.64 

0.9214 

0.3 

0.0251 

Adhes- 
ive 

Test 
specimen 

0.5 

0.5 

0.192 

0.3 

0.3 

Table  4.4  Geometric  parameters  of  baseline  design  and  test  specimen. 


Panel 

Design 

Test  specimen 

Panel  length  (in) 

30 

32 

Panel  width  (in) 

32 

32 

Blade  height  (in) 

3.0705 

3.0723 

Skin  ply  thickness  (in) 

0.00520 

0.00555 

Blade  ply  thickness  (in) 

0.00520 

0.00503 

Flange  ply  thickness  (in) 

0.00520 

0.00491 

Adhesive  thickness  (in) 

0.0 

0.0121 
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Table  4.5  Buckling  load  factors  obtained  from  STAGS  with  various  boundary  conditions 

and  material  properties. 


Boundary  conditions 

shell  unit 

shell  unit 

\^  ^  1  ClCIIIClll^ 

4  edges  simply  supported 
(baseline  design) 

1.166 

1.296 

2  edges  free,  2  edges  clamped 
(baseline  design) 

1.197 

1.340 

2  edges  free,  2  edges  clamped 
(potted  region,  test  material,  32-inch 
panel  length) 

1.154 

1.286 

*  Incluc 

es  shear  deformation. 

Figure  4.4  The  STAGS  predicted  buckling  mode  shape  of  the  perfect  baseline  design 
panel  subjected  to  uniaxial  compression  (load  factor=1.166). 
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Figure  4.5  The  predicted  buckling  mode  shape  of  the  blade  stiffened  panel  with  9-node 
shell  elements  from  the  ABAQUS  (load  factor=1.218). 


Figure  4.6  The  predicted  buckling  mode  shape  of  the  blade  stiffened  panel  from  the 

STAGS  (load  factor=  1.154). 
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An  examination  of  tiie  STAGS  model  in  Figure  4.3  shows  that  there  is  double- 
counting  of  material  between  the  blade  stiffener  and  the  flange-skin  combination.  This  is 
because  of  the  way  the  nodes  are  located  in  the  blade  elements  and  in  the  element  that 
represents  the  flange-skin  combination.  Both  elements  have  mid-plane  nodes  leading  to 
an  overlap  equal  to  the  thickness  of  the  blade  with  a  width  equal  to  half-thickness  of  the 
flange-skin  combination.  This  overlap  is  avoided  in  the  modeling  approach  described 
earlier  for  the  ABAQUS  model.  This  additional  material  due  to  the  overlap  in  the  model 
is  expected  to  increase  the  pre-buckling  axial  stiffness  of  the  panel.  In  order  to  estimate 
this  increase  the  STAGS  model  was  modified  similarly  to  the  ABAQUS  model,  and  the 
linear  buckling  analysis  was  performed  on  the  modified  model.  The  results  are  presented 
in  Table  4.6  and  are  also  compared  with  the  experimental  results.  The  increase  in  area 
due  to  the  overlap  contributed  to  about  one  half  of  the  total  stiffness  difference  between 
the  analysis  and  test  results. 


Table  4.6  Comparison  of  prebuckling  stiffness  and  buckling  load  factors 

(STAGS  results). 


EA/L  (kip/in) 

EA  (kip) 

Buckling  load  factor 

Model  with 
overlap 

3931.2 

125,800 

1.15 

Model  with 
no  overlap 

3684.4 

117,900 

1.23 

Experiment 

3453.9 

110,500 

1.09 
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4.5.3  Summary  of  Differences  between  Design  Model  and  Test  Model 

In  summary,  the  PASCO  model  used  for  designing  the  panel  had  several  modeling 
simplifications  and  compensating  factors;  their  effects  are  listed  in  Table  4.7.  The  two 
major  model  simplifications  were: 

1.  PASCO  does  not  account  for  shear  deformation,  which  is  significant  for  thick 
composite  panels,  and  this  reduces  the  buckling  load  by  11%. 

2.  PASCO  employs  simple  support  boundary  conditions.  The  difference  due  to 
boundary  conditions  was  only  about  1%  because  the  buckling  mode  is  local. 

To  obtain  a  more  robust  design,  the  PASCO  model  was  subjected  to  an  additional 
shear  load  and  3%  imperfection.  Both  had  substantial  effects  on  the  design  load.  However, 
because  the  panel  was  designed  only  for  imperfection  of  one  sign,  it  became  more  sensitive 
to  an  imperfection  of  the  opposite  sign.  Finally,  because  of  the  substantial  thickness  of  the 
blade,  it  was  also  found  that  a  centerline  modeling,  which  is  common  in  thin  walled 
structures,  produces  about  seven  percent  increase  in  the  prebuckling  stiffness,  and  a  seven 
percent  reduction  in  buckling  load.  The  opposite  effects  are  explained  by  the  fact  that  the 
overlap  draws  more  loads  into  the  blade,  which  is  the  critical  element.  Overall,  the  more 
accurate  analytical  model  predicts  a  buckling  load  higher  by  about  18%  than  the  design 
load;  however,  this  does  not  take  into  account  any  imperfections.  The  actual  buckling  load 
was  only  about  10%  higher  than  the  design  load.  It  can  be  concluded  that  for  this  panel  the 
simplified  model  used  in  PASCO,  together  with  the  shear  and  imperfection  loading  added 
for  robustness,  worked  reasonably  well.  The  other  two  designs  that  were  tested  also  buckled 
at  or  slightly  above  the  design  load 
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Table  4.7  Differences  between  baseline  design  analysis  and  actual  test  panel  analysis. 


Baseline  design 

Test  panel 

Effect  on  buckling 
load 

Loading  (lbs/inch) 
Axial  compression 
Inplane  shear  load 

20,000 
5,000 

20,000 
0 

-20% 

Imperfection 
(initial  bow  type) 

-3%  of  panel  length 

Not  measured 

-30%  (with  shear 
load) 

-18%  (w/o  shear  load) 

Boundary  Condition 
Loaded  edges 
Unloaded  edges 

Simple  support 
Simple  support 

Clamped  (potted) 
Unsupported 

1%  (with  test 
material) 

Transverse  shear 
deformation 

no 

yes 

11% 

4.6  Nonlinear  Analysis 

Although  the  linear  buckling  loads  provide  a  measure  of  the  compressive  load 
carrying  capacity  of  the  stiffened  panels,  the  test  results  indicate  that  the  panels 
underwent  substantial  nonlinear  transverse  deformations  prior  to  failure.  Hence  it  was 
decided  to  perform  a  nonlinear  analysis  in  order  to  understand  the  effects  of  boundary 
conditions  including  the  eccentricity  in  load  application.  The  nonlinear  analysis  was 
started  without  applying  any  initial  imperfection,  but  the  differences  in  the  stacking 
sequence  and  the  differences  in  material  properties  between  the  skin  and  blades  induce 
bending  deformation.  The  modified  Riks  path  following  algorithm  in  STAGS  was  used 
for  the  nonlinear  analysis.  The  computation  time  required  for  the  nonlinear  analysis  was 
an  order  of  magnitude  higher  than  for  the  linear  analysis,  indicating  significant 
nonlinearity  (probably  near  the  buckling  load).  In  the  following  sub-sections  we  discuss 
the  results  of  the  nonlinear  analyses,  and  make  a  comparison  between  experimental  and 
analytical  results. 
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4.6. 1  Compressive  Load  Versus  End-shortening 

The  compressive  load  versus  end-shortening  deflection  curves  from  the 
STAGS  nonlinear  static  analyses,  the  test,  and  a  linear  fit  (regression  analysis)  of  the 
measured  data  are  shown  in  Fig.  4.7.  Recall  that  the  test  panel  designated  as  GA2461 
in  Reference  [80]  is  the  baseline  design  for  the  present  paper.  In  comparison  to  the 
baseline  panel,  two  other  test  panels  (GA2414,  GA2458)  from  Reference  [80]  have 
slightly  different  geometries  and  stacking  sequences.  As  expected,  their  compressive 
load  versus  end-shortening  curves  from  the  tests  exhibit  a  similar  trend  except  at  the 
initial  stage  of  loading.  The  prebuckling  stiffness  (EA)  was  calculated  from  the  slopes 
of  the  linear  portions  of  the  experimental  load  versus  end-shortening  curves  for  the 
three  panels  and  by  multiplying  the  slopes  by  the  panel  length.  The  prebuckling 
stiffness  from  the  experiments  is  compared  in  Table  4.8  with  the  prebuckling  stiffness 
predicted  using  STAGS  for  test  panel  GA2461.  It  is  observed  that  the  prebuckling 
stiffness  of  the  test  panels  is  about  6%  lower  than  the  analytical  value. 


Table  4.8  Comparison  of  the  prebuckling  stiffness  and  buckling  load  results  from 


STAGS  and  experiments. 


EA/L  (kip/in) 

EA  (kip) 

Buckling  load 
factor 

STAGS  nonlinear 
(GA2461) 

3931.2 

117,900 

1.23 

Experimental 
Result  (GA261) 

3453.9 

110,500 

1.09 

Experimental 
Result  (GA2414) 

3347.8 

107,100 

0.94 

Experimental 
Result  (GA2458) 

3333.8 

106,700 

1.07 
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Figure  4.7  Compressive  load  versus  end-shortening  from  analysis  and  experiments. 


Load  vs.  Out-of-Plane  displacement 
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Figure  4.8  Load  versus  out-of-plane  displacements  of  the  skin  at  the  center  bay. 
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Figure  4.9  Load  versus  out-of-plane  displacements  of  the  stiffener  at  selected  location  in 

Fig.2. 


Figure  4.10  Variation  of  the  out-of-plane  displacements  along  the  panel  lengthwise 

direction. 
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4.6.2  Compressive  Load  Versus  Out-of  -Plane  Deformations 

The  layout  of  the  DCDTs  used  to  measure  displacements  in  the  test  panel  is 
shown  in  Figure  4.2.  The  out-of-plane  displacements,  measured  from  DCDTs  at  selected 
locations  in  Figure  4.2,  are  shown  in  Figures.  4.8  and  4.9.  The  results  in  Figure  4.8 
show  that  out-of-plane  displacements  were  initiated  at  an  early  stage  of  the  loading  and 
increase  linearly  in  proportion  to  the  loading.  This  observation  suggested  the  possibility 
of  loading  eccentricities  along  the  load  introduction  edge  or  rigid  body  rotation  of  the 
panel  with  respect  to  the  clamped  edge  in  addition  to  the  effects  of  geometric 
imperfections.  Figure  4.9  shows  that  the  out-of-plane  displacements  of  the  blade 
stiffeners  were  also  started  at  an  early  stage  of  the  loading.  Except  DCDT  11,  the  out-of- 
plane  deformations  were  an  order-of-magnitude  lower  than  those  of  the  skin  in  Figure 
4.8.  Furthermore,  significant  nonlinear  response  was  only  exhibited  near  the  failure  load. 
The  large  nonlinear  response  of  DCDTl  1  throughout  the  axial  loading  was  probably  due 
to  the  effect  of  the  unsupported  side-edge  boundaries.  The  out-of-plane  displacement 
variations  along  the  length  of  the  panel  (DCDTs  1-4  in  Figure  4.2)  at  selected  load  levels 
are  shown  in  Figure  4.10.  The  results  in  Figure  4.10  indicate  that  bending  occurred  in 
the  test  panel  in  addition  to  the  end  shortening  due  to  the  axial  compressive  load.  The 
load  versus  out-of-plane  displacements  across  the  panel  mid-length  can  be  found  in  Ref. 
[83]. 

In  order  to  explain  the  substantial  prebuckling  bending,  a  combination  of 
different  geometric  imperfections  and  loads  applied  at  a  small  angle  to  the  axial 
direction  were  analyzed  to  determine  their  influence  on  the  observed  out-of-plane 
displacements.  The  capability  of  STAGS  to  model  geometric  imperfections  in  the  shape 
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of  the  buckling  modes  was  used  for  these  analyses.  Various  combinations  of 
imperfection  amplitudes  and  load  angles  were  considered.  Although  for  some 
combinations  we  could  reproduce  the  test  results  partially  [83],  obtaining  the  right 
imperfection  and  the  load  introduction  angle  seemed  elusive.  This  difficulty  suggests 
that  we  look  elsewhere  for  the  source  of  the  prebuckling  bending. 

4.6.3  Contact  Between  the  Panel  and  Loading  Platen 

Hilburger  [97]  investigated  the  effects  of  non-uniform  load  introduction  and 
boundary  condition  imperfections  on  the  compression  response  of  composite 
cylindrical  shells  with  cutouts.  He  defined  the  non-uniform  load  distribution  as 
anything  other  than  uniform  axial  displacement  of  cylinder's  loading  surface  and 
found  two  sources  of  non-uniform  load  introduction.  One  was  due  to  lack  of  planarity 
in  the  loading  surfaces  of  the  specimen  and  the  loading  platens.  The  other  source  was 
due  to  tilt  of  the  loading  platen  with  respect  to  the  specimen  before  the  loading  began. 
He  measured  the  top  and  bottom  loading  surface  imperfections  as  well  as  potting 
thickness.  Then  the  imperfection  data  was  fit  to  curves  and  input  into  the  STAGS 
models.  Furthermore,  the  test  frame  loading  platen  was  modeled  as  rigid  flat  plates  and 
generalized  contact  definitions'  given  in  STAGS  were  used. 

A  similar  modeling  approach  was  used  in  the  present  study  in  order  to  identify 
the  causes  of  the  substantial  out-of-plane  deformations  in  addition  to  nonlinear  end 
shortening  during  the  early  stage  of  the  test.  The  loading  platen  was  modeled  as  a  rigid 


•  Generalized  contact  definition  means  that  contact  points  are  calculated  by  STAGS  rather  than  specified  by  the  user. 
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flat  plate  in  the  STAGS  analysis.  Because  the  loading  surface  imperfections  were  not 
measured  before  the  test,  they  were  not  considered  in  this  study.  Instead,  it  was  assumed 
that  the  rigid  loading  platen  was  initially  contacted  at  the  tip  of  the  blade.  Stiffeners  were 
assumed  to  have  a  small  tilt  angle  with  respect  to  the  load  introduction  edge  of  the  test 
specimen,  as  shown  in  Figure  4. 1 1 . 

STAGS  uses  generalized  contact  definitions  to  check  for  contact  and  to  construct 
actual  contact  elements  coupling  contacting  points  with  contacted  shell  elements  as  the 
analysis  progresses.  In  doing  this,  STAGS  uses  penalty  functions  to  enforce  a 
displacement-compatibility  constraint  between  each  contact  point  and  each  element  with 
which  it  is  in  contact.  STAGS  utilizes  analyst-supplied  stiffness  versus  displacement 
information  to  compute  the  forces  resulting  from  the  small  contact-surface  penetration 
that  may  occur.  A  contact  element  is  conceptually  a  nonlinear  spring  connecting  the 
contact  point  to  the  surface  of  the  contacted  element.  This  nonlinear  spring  typically  has 
a  low  stiffness  and  generates  a  small  force  when  the  contact-surface  penetration  is  small, 
but  it  gets  progressively  stiffer  and  generates  a  larger  force  as  the  penetration  increases 
[23]. 

Generalized  contact  definitions  were  implemented  to  the  finite  elements  that 
simulate  the  rigid  platen  and  the  load  introduction  edge  of  the  test  specimen.  The 
selection  of  proper  stiffnesses  of  the  contact  elements  for  present  analyses  is  rather 
arbitrary.  Thus,  several  nonlinear  analyses  were  performed  to  simulate  the  observed  out- 
of-plane  deformation  response  of  the  test  specimen  by  changing  both  tilt  angles  and  the 
stiffnesses  of  contact  elements  between  the  loading  platen  and  load  introduction  edge. 
The  combinations  of  the  tilt  angles  and  stiffnesses  of  the  contact  elements  used  for  the 
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analyses  are  summarized  in  Table  4.9,  and  the  corresponding  load  versus  end 
shortening  results  are  shown  in  Figure  4.12.  The  results  in  Figure  4.12  suggest  that  the 
computed  end-shortening  response  strongly  depend  on  the  user-supplied  input  data  in 
Table  4.9.  The  response  of  Model  9  was  the  closest  to  that  of  the  test  panel  in  Figure 
4.12.  The  computed  out-of-plane  displacements  of  the  skin  and  stiffeners  of  Model  9  at 
the  selected  DCDTs  location  are  further  examined  as  shown  in  Figures.  4.13  and  4.14, 
respectively.  It  is  observed  that  the  response  of  Model  9  in  Figure  4.13  shows  good 
correlation  with  those  of  the  skin  of  the  test  panel  in  Figure  4.8  during  the  early  stages  of 
the  load.  However,  Model  9  exhibits  considerable  nonlinear  behavior  of  the  out-of- 
plane  deformations  of  skin  when  the  compressive  load  is  above  400,000  lb,  which  was 
not  present  in  the  response  of  the  test  panel  in  Figure  4.8. 


Table  4.9  Summary  of  the  tilt  angles  and  the  stiffness  of  contact  element. 


Panel  length 
(inch) 

Tilt  angle 
(degree) 

Stiffness  of  contact  element 
(lb/inch) 

Model  7 

30 

0.01 

Disp. 

0.005 

0.05 

0.1 

0.2 

1.0 

Force 

1.0e3 

1.0e5 

l.OeV 

2.0e8 

2.0e8 

Model  8 

30 

0.01 

Disp. 

0.0001 

0.03 

0.05 

0.2 

1.0 

Force 

1.0e3 

l.Oe? 

l.OeS 

2.0e8 

2.0e8 

Model  9 

30 

0.01 

Disp. 

0.005 

0.03 

0.05 

0.2 

1.0 

Force 

1.0e3 

1.0e5 

l.OeS 

2.0e8 

2.0e8 

Model  10 

30 

0.005 

Disp. 

0.005 

0.03 

0.05 

0.2 

1.0 

Force 

1.0e3 

1.0e5 

1.0e8 

2.0e8 

2.0e8 

The  computed  load  versus  out-of-plane  deformation  of  the  stiffener  at  location 
of  the  DCDT  11  shows  significant  nonlinear  behavior  in  Figure  4.14.  This  significant 
nonlinear  response  was  also  observed  from  the  measured  response  of  the  actual  test 
panel  in  Figure  4.9.  In  general,  the  finite  element  model  with  the  generalized  contact 
definitions  improved  correlation  between  the  measured  and  predicted  out-of-plane 
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deformation.  However,  the  details  of  the  displacements  are  considerably  different.  It  is 
concluded  that  some  combination  of  tilt  angles  and  the  contact  stiffnesses  can  produce 
the  observed  pattern,  but  there  may  be  some  other  contribution  to  the  out-of-plane 
displacements. 


I 

i 


TILT 
ANGLE 


Figure  4. 11  Schematic  of  blade  stiffened  panel  and  loading  platen. 
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Figure  4.12  Compressive  load  versus  end-shortening  from  analysis  with  contact  models 
and  experiment.  Model  numbers  refer  to  Table  9. 
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Figure  4.13  Load  versus  out-of-plane  displacements  of  the  skin  at  the  center  bay  from 

Model  9  analysis. 
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Figure  4.14  Load  versus  out-of  plane  displacements  of  the  stiffeners  at  selected  locations 

from  Model  9  analysis. 


4.7  Conclusion 
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Analytical  models  using  several  structural  analysis  models  were  used  to  assess 
the  adequacy  of  the  design  model  and  the  correlation  with  experimental  results  for  a 
stiffened  panel  designed  using  the  PASCO  program.  Of  the  effects  neglected  by  the 
simple  model,  shear  deformation  was  the  most  important,  accounting  for  about  11% 
difference  in  buckling  load.  The  effect  of  simplified  (simple  support)  boundary 
conditions  was  small.  The  addition  to  the  design  model  of  shear  loads  and  imperfections 
to  improve  the  robustness  of  the  result  did  help,  even  though  the  inclusion  of  one-sided 
imperfection  apparently  induced  sensitivity  to  imperfection  of  the  opposite  sign.  Overall, 
the  simplified  model  did  produce  designs  that  in  the  experiments  failed  slightly  above 
the  design  load. 

The  most  significant  difference  between  the  analytical  predictions  and 
experimental  measurements  was  the  substantial  out-of-plane  pre-buckling  deformations. 
To  explain  these  differences,  imperfections,  load  eccentricities,  and  loading  platen  tilt 
angles  were  considered.  Of  these  the  loading  platen  tilt  produced  similar  patterns  of 
deformation,  but  these  had  more  nonlinear  characteristics  than  the  measured 
deformations. 


CHAPTER  5 

BUCKLING  AND  POSTBUCKLING  ANALYSIS  OF  A  STIFFENED  PANEL  WITH 

SKIN-STIFFENER  DEBOND 

5.1  Introduction 

In  contrast  to  most  of  simplified  approaches  discussed  in  literature  survey  in 
Chapter  1 ,  finite  element  based  approach  can  be  applicable  with  high  fidelity  to  more 
general  and  realistic  structures.  Therefore,  the  axial  compressive  behavior  of  a  stiffened 
panel  with  skin-stiffener  debond  was  explored  in  this  chapter  by  employing  finite 
element  model.  Two  different  finite  element  models,  where  nodes  of  the  panel  skin 
elements  and  the  stiffener  flange  elements  are  located  on  the  mid-plane  of  each  element 
or  interface  between  the  skin  elements  and  flange  elements  with  offset,  were  used.  The 
nodes  corresponding  to  the  top  of  the  skin  and  bottom  of  flange  are  connected  with 
either  the  elastic  spring  fastener  elements  or  multi-point  constraint  equations. 

In  order  to  verify  present  finite  element  modeling  approach,  laminated 
composite  plates  with/without  the  through-the-width  delamination  were  first  modeled. 
Both  single  delamination  and  multiple  delaminations  were  considered.  Then,  stiffened 
composite  panels  with  skin-stiffener  debond  were  modeled.  Stiffened  composite  panels 
with  a  single  stiffener  as  well  as  multiple  stiffeners  were  also  considered.  Buckling  and 
postbuckling  analyses  were  conducted  using  STAGS.  Comparison  was  made  with 
available  buckling  analysis  results.  Next,  numerical  examples  of  computing  energy 
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release  rate  in  the  context  of  plate  as  discussed  in  Chapter  3  are  given  for  predicting 
debond  extension  using  the  strain  energy  derivative  method,  the  virtual  crack  closure 
technique,  and  the  crack-tip  force  method. 

5.2  Finite  Element  Model 

Two  finite  element-modeling  approaches  for  the  stiffened  panel  are  commonly 
used  in  the  literature.  One  approach  is  to  model  the  skin  with  plate  elements  and  to 
model  the  stiffener  with  beam  elements.  The  other  is  model  both  skin  and  stiffener  with 
plate  elements.  The  second  modeling  approach  appears  more  attractive  for  modeling 
debonded  region  between  skin  and  flange,  and  it  was  chosen  in  this  study.  Furthermore, 
two  different  finite  element  models,  designated  as  Model  I  and  Model  II,  respectively, 
were  also  considered  (see  Figure  5.1).  Both  skin  and  flange  elements  of  Model  I  have 
offset  nodes  with  small  gap  at  the  interface  region,  while  nodes  in  Model  II  are  located 
on  the  mid-planes  of  skin  and  flange,  respectively.  In  Model  I,  in  order  to  satisfy 
compatibility  conditions  of  intact  interface  nodes  located  directly  above  the  skin  and 
below  the  flange,  each  nodal  degree  of  freedom  was  constrained  with  elastic  spring 
fastener  elements  with  very  high  spring  constants.  Nodes  located  on  debonded 
interface  between  skin  and  flange  were  connected  with  elastic  fastener  elements,  which 
have  only  an  axial  degree  of  freedom  with  very  high  stiffness  in  compression  and  zero 
stiffness  in  tension.  This  can  prevent  physically  unrealistic  nodal  penetration  between 
skin  and  flange  during  the  postbuckling  analysis.  Friction  against  sliding  of  the 
debonded  surface  was  not  considered.  In  Model  II,  Multi-point  constraints  were 
imposed  at  the  interface  of  intact  skin  and  corresponding  flange  nodes  to  satisfy  the 
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displacement  compatibility  condition  as 


u 


-  V 


/ 


(5.1) 


where  u,  v,  and  w  are  displacements  in  x,  y,  and  z  direction,  %  and     are  rotations,  h 
is  thickness,  and  subscripts  s  and/  denote  skin  and  flange,  respectively.  In  the 
debonded  region,  the  same  modeling  approach  as  Model  I  was  used. 

One  advantage  of  model  I  over  model  II  is  such  that  displacement 
compatibility  conditions  at  interface  nodes  in  Equation  (5.1)  is  not  necessary.  The  other 
advantage  of  model  I  is  that  computation  of  displacements  behind  crack  tip  is  easy. 
However,  Model  I  cannot  handle  multiple  delaminations  located  through-the-thickness 
direction  while  Model  II  can. 


Figure  5.1  Two  different  modeling  approach  for  blade  stiffened  panel. 
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5.3  Buckling  Analysis  Results  of  Multiple  Delaminated  Plate 


A  plane  woven  fabric  glass  fiber  reinforced  composite  panel  with  three 
through-the-width  delaminations  located  in  the  middle  of  the  plate  from  Suemasu  [56] 
(as  shown  in  Figure  5.2)  was  examined  using  Model  II  approach.  The  material 
properties  are  shown  in  Table  5.1. 

The  STAGS  finite  element  model  for  the  plate  has  a  total  of  4  branched  shell 
units  with  four  node  shell  elements  (element  410),  and  each  branched  shell  unit  has  41 
by  11  nodes.  Total  of  297  elastic  fastener  elements  (element  130)  was  used  to  prevent 
penetration  of  contact  surfaces  in  debonded  region  of  the  plate.  The  compressive 
stiffness  of  an  elastic  fastener  element  used  in  this  study  is  113  MN/m.  In  the  intact 
region  of  plate,  a  total  of  3168  constraint  equations  was  used  to  satisfy  the 
compatibility  conditions  of  shell  unit  interfaces.  Axial  compressive  load  (2500  N)  was 
applied  to  load  introduction  edge  with  uniform  end-shortening  constraint.  In  order  to 
ensure  a  uniform  stress  condition  for  the  entire  panel  length,  an  incremental  boundary 
condition  option  was  chosen  to  prevent  axial  bending  during  the  prebuckling  stage. 
The  computed  buckling  load  and  buckling  load  from  Reference  [56]  are  given  in  Table 
5.2.  The  first  and  second  buckling  mode  shapes  from  Reference  [56]  and  present 
analysis  are  shown  in  Fig.  5.3  and  Fig.  5.4,  respectively.  It  is  seen  that  present  buckling 
analysis  results  agree  well  with  the  corresponding  buckling  analysis  results  from 
Reference  [56].  However,  the  buckling  load  from  the  experiment  is  lower  than  the 
computed  buckling  load.  Suemasu  [56]  indicated  that  insufficient  clamped  condition 
during  his  experiment  might  be  responsible  for  low  buckling  load. 
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Table  5.1  Plane  woven  fabric  fiber  reinforced  laminates  material  properties  [Ref.  56]. 


Young's  modulus  (longitudinal) 

El  =  20.2  GPa 

Young's  modulus  (transverse) 

E2  =  21.0  GPa 

Young's  modulus  (through-thickness) 

E3=  10.0  GPa 

Shear  modulus  (inplane) 

Gi2  =  4.15  GPa 

Shear  modulus  (through-thickness) 

Gi3  =  4.0GPa 

Poisson's  ratio  (inplane) 

Vi2  =  0.16 

Poisson's  ratio  (through-thickness) 

Vi3  =  0.3 

Figure  5.2  Glass  fiber  reinforced  composite  specimen  [Ref.  56]. 
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Table  5.2  Comparison  of  buckling  load  for  delaminated  plates  (  Pcr=2500  N). 


Buckling  load  factor 

Present 

Suemasu  [56] 

Mode  I 

0.440 

0.432 

Mode  II 

0.801 

0.795 

Figure  5.3  Buckling  modes  of  the  plate  with  three  delamaninations  from  Ref.  [56]. 


Figure  5.4  Computed  buckling  modes  using  model  II  approach, 
(a)  First  lowest  buckling  mode,  (b)  Second  lowest  buckling  mode. 
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5.4  Buckling  and  Postbuckling  Analysis  Results  of  Plate  with  Single  Delamination 

In  order  to  see  the  effects  of  the  delamination  location  through  thickness,  a 
unidirectional  graphite/epoxy  laminated  plate  with  single  through-the-width 
delamination  was  considered  (see  Figure  5.5).  The  geometry  and  material  properties 
used  for  the  graphite/epoxy  laminates  are  given  in  Table  5.3.  Axial  compressive  load 
(1,250  lb)  was  applied  to  load  introduction  edges.  Three  different  thickness  ratios  of  top 
and  bottom  subraminates  were  considered  and  buckling  load  factors  are  summarized  in 
Table  5.4.  The  computed  buckling  modes  agree  well  with  deformed  shapes  obtained 
from  the  test  in  Reference  [53]  (see  Figures  5.6,  5.7).  As  expected,  the  delamination 
located  near  the  surface  of  the  plate  has  the  lowest  buckling  load,  where  local  buckling 
is  the  dominant  buckling  mode.  As  the  thickness  ratio  of  sublaminates  increases,  the 
buckling  mode  changes  from  local  buckling  to  mixed  buckling  mode  (see  Figure  5.8) 
and  then  from  mixed  to  global  buckling  mode. 


Thickness  of  top  sublaminate  (hi) 
Thickness  of  bottom  sublaminate  (h2) 


Figure  5.5  Graphite/epoxy  laminated  plate  with  single  through-the-width 
delamination.  All  dimensions  are  inch. 
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TABLE  5.3  Graphite  epoxy  lamina  material  properties 


Young  s  modulus  (longitudmal) 

El  =19.50  X  lO^psi 

Young's  modulus  (transverse) 

E2=1.48  X  lO^psi 

Shear  modulus 

Gi2=0.80x  lO^psi 

Poisson's  ratio 

Vi2  =  0.3 

Plate  length  x  width 

5.0  X  1.25  in. 

Plate  total  thickness 

0.15  in. 

Delamination  length 

2.5  in. 

Table  5.4  Buckling  load  delaminated  plates  with  different  thickness  ratios  of 


sublaminates  (Pappl.=] 

1,250  lb). 

hl/h2=0.5 

hl/h2=0.8 

hl/h2=1.0 

Buckling  load  factor 

3.44 

5.42 

5.79 

(a) 


Figure  5.6  Global  buckling  mode  from  test  and  analysis, 
(a)  Buckling  mode  in  mid-plane  delamination  [53]. 
(b)  Buckling  mode  from  analysis  (hl/h2=1.0). 


(a) 


1  u-ai 


XL. 


7  I 


(b) 


Figure  5.7  Local  buckling  mode  from  test  and  analysis, 
(a)  Buckling  mode  in  near-surface  delamination  [53]. 
(b)  Buckling  mode  from  analysis  (hl/h2=0.5). 
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Figure  5.8  Mixed  buckling  mode  from  analysis  (hl/h2=0.8). 


Figure  5.9  Load  versus  end-shortening  with  various  thickness  ratios  of  sublaminates. 
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Postbuckling  analysis  was  also  conducted  to  see  the  load  carrying  capacity  of  the 
delaminated  plates.  Imperfection  was  assumed  in  the  form  of  the  first  mode  shape 
obtained  from  linear  buckling  analysis  and  the  maximum  magnitude  of  imperfection  for 
postbuckling  analysis  was  1  %  of  the  top  sublaminte  thickness.  The  modified  Riks 
algorithm  in  STAGS  was  used  for  this  computation.  Loads  versus  end-shortening  curves 
with  different  thickness  ratios  of  sublaminates  are  shown  in  Figure  5.9.  From 
postbuckling  analysis  results  of  Figure  5.9,  one  can  find  that  the  ultimate  load  carrying 
capacity  of  the  composite  plate  with  a  small  thickness  ratio  of  sublaminates  (e.g., 
hl/h2=0.5)  is  significantly  higher  than  the  buckling  load  obtained  from  linear  bifurcation 
analysis.  This  trend  was  also  observed  from  the  test  results  in  Reference  [53]. 

5.5  Buckling  Analvsis  Results  Stiffened  Panel  with  Debond 

Two  blade  stiffened  panels,  available  from  References  [83]  and  [86], 
respectively,  were  examined  in  detail.  The  configurations  of  blade  stiffened  panel, 
material  properties,  and  stacking  sequences  from  References  [83]  and  [86]  are 
summarized  in  Table  5.5.  The  computed  energy  release  rate  data  is  available  but  the 
buckling  load  is  not  available  in  Reference  [86].  Conversely,  experimentally  determined 
buckling  load  data  for  panel  without  skin-stiffener  debond  is  available  but  energy  release 
rate  is  not  available  in  Reference  [83].  Thus,  both  panels  with  single  blade  stiffener  were 
considered.  The  panel  with  a  single  blade  stiffener  (Figure  5.1 1)  has  a  total  of  544  nine- 
node  shell  element  (Element  480)  with  2341  nodes,  two  node  fastener  elements 
(Element  130)  were  used  for  the  joining  skin  and  corresponding  flange  nodes.  An  axial 
compressive  load  of  3.503MN/m  (20,0001b/in)  for  the  panel  from  Reference  [83]  was 


TABLE  5.5  Graphite  epoxy  lamina  material  properties  from  Ref.  [83]  and  Ref.  [86] 


Ref.  [86] 

Ref.  [83] 

Young's  modulus 
(longitudinal) 

E,  =19.50  X  10^  psi 

E,  =18.50  X  10* 
psi 

Young' s  modulus 
(transverse) 

E2=1.48  X  lO^psi 

E2=1.48  X  10*  psi 

Shear  modulus 

G,2  =0.80  X  10^  psi 

G,2  =0.80  X  10* 
psi 

Poisson's  ratio 

Vi2  =  0.3 

v,2  =  0.3 

TABLE  5.6  Geometric  parameters  of  a  blade  stiffened  panel  from  Ref.  [83]  and 
Ref.  [86] 


Panel 

Ref.  [86] 

Ref.  [83] 

Panel  length  (in) 

21 

30 

Panel  width  (in) 

5 

8 

Blade  height  (in) 

1.45 

3.0705 

Skin  thickness  (in) 

0.083 

0.208 

Blade  thickness  (in) 

0.105 

0.3536 

Flange  width  (in) 

2.0 

2.4 

Flange  thickness  (in) 

0.067 

0.3536 
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applied  with  uniform  end-shortening  constraint.  In  load  introduction  edges,  Poisson 
expansion  was  allowed.  Both  symmetric  and  free  boundary  conditions  of  unloaded  side 
edges  were  examined.  The  symmetric  boundary  conditions  represent  infinitely  wide 
panel  but  with  debonded  stiffeners  at  uniform  spacing.  Rigid  body  motion  was 
constrained  along  one  side  edge. 

5.5.1  Buckling  analysis  results  of  debonded  stiffened  panels 

Buckling  mode  of  perfect  panel  with  symmetry  boundary  condition  agree  well 
with  buckling  mode  obtained  in  Reference  [83]  (see  Figures  5.10  and  5.11).  Figures 
5.12-5.13  show  buckling  modes  of  debonded  stiffened  panels.  As  the  length  of  the 
debond  increases,  the  buckling  mode  changes  from  global  buckling  mode,  where  blade 
tip  buckling  is  dominant,  to  mixed  buckling  mode  and  local  buckling  mode, 
respectively.  This  mode  transition  is  expected  because  the  axial  and  bending  stiffness  of 
the  blade-flange  combination  is  higher  than  the  stiffness  of  skin.  Thus,  the  buckling 
mode  of  the  debonded  skin  is  similar  to  that  of  one-dimensional  thin-film  analysis.  This 
suggests  that  one  dimensional  beam-plate  model  may  be  still  useful  to  predict  buckling 
load  of  this  problem  when  local  buckling  is  dominant. 

In  order  to  see  the  effect  of  free  side  edges  on  the  buckling  load,  the  two 
aforementioned  cases  of  boundary  conditions  were  examined,  and  buckling  load 
variations  with  respect  to  debond  ratio  (debond  length  divided  by  panel  length)  are 
shown  in  Figure  5.14.  Based  on  the  results  in  Figure  5.14,  it  can  be  seen  that  there  is  a 
critical  debond  length  which  does  not  reduce  noticeably  the  buckling  load,  and  this 
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critical  debond  length  varies  with  the  boundary  conditions.  Furthermore,  local  buckling 
is  relatively  insensitive  to  boundary  conditions  compared  to  global  buckling.  Variations 
of  buckling  loads  with  respect  to  different  stiffener  geometry  and  stacking  sequence  for 
the  panel  configuration  from  Ref.  86  are  shown  in  Figures  5.15  and  5.16,  respectively. 
As  expected,  buckling  loads  increase  in  proportion  to  the  increase  of  stiffener  thickness. 
When  the  debond  ratio  is  greater  than  about  0.4,  contributions  of  increased  stiffener 
thickness  on  buckling  load  are  relatively  small  due  to  the  local  buckling.  Figure  5.16 
shows  that  the  panel  with  all  0  degree  plies  shows  small  buckling  loads  compared  to  the 
panels  with  cross  plies  or  balanced  symmetry  stacking  sequences  when  debond  ratios  are 
less  than  0.2.  This  suggests  that  stiffened  panel  with  all  0  degree  plies  can  lose  stability 
with  ease  in  spite  of  high  axial  stiffness.  Buckling  analysis  of  the  corresponding  panel 
without  debond  was  also  conducted  using  PANDA2  [15-17].  The  predicted  buckling 
load  factors  from  PANDA2  with  and  without  neglecting  the  redistribution  of  axial 
compressive  load  due  to  local  buckling  are  0.81  and  0.92,  respectively.  Those  buckling 
load  factors  are  fairly  close  to  finite  element  analysis  results,  0.95,  in  Figure  5.16. 


Figure  5.10  Buckling  mode  of  stiffened  panel  from  Reference  [83] 


Figure  5.1 1  Global  buckling  mode. 


Figure  5. 12  Mixed  buckling  mode. 


Figure  5.13  Local  buckling  mode. 
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Figure  5.14  Variations  of  buckling  load  with  respect  to  debond  ratio. 
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Figure  5.15  Effects  of  the  flange  geometry  on  buckling  loads 


Figure  5.16  Effects  of  the  Stacking  Sequences  on  Buckling  Loads 
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5.5.2  Postbuckling  analysis  results 

As  mentioned  previously,  stiffened  laminated  composite  flat  panels  usually 
exhibit  stable  postbuckling  behavior,  which  may  lead  to  significant  differences  between 
buckling  load  and  ultimate  failure  load.  If  the  structure  exhibits  considerable  nonlinear 
prebuckling  behavior  due  to  initial  imperfection  or  excessive  bowing  associated  with 
local  buckling,  then  it  is  necessary  to  perform  nonlinear  analysis. 

A  nonlinear  analysis  was  performed  for  the  panel  made  up  of  all  0  degree  ply 
laminates  from  Reference  [86].  An  axial  compressive  load,  22,24 IN  (50001b),  was 
applied  with  uniform  end-shortening  constraint.  Initial  imperfection  was  assumed  as  first 
mode  shape  obtained  from  linear  buckling  analysis.  Load  versus  end-shortening  curves 
with  various  debond  ratios  are  shown  in  Figure  5.17.  Figures  5.18-19  show  loads 
versus  out-of-plane  deformations  with  several  magnitudes  of  imperfection  and  debond 
ratios,  respectively.  The  load  versus  out-of-plane  deformation  exhibits  stable  nonlinear 
postbuckling  behavior.  This  suggests  that  linear  buckling  analysis  results  of  these 
specific  problems  can  provide  overly  conservative  estimation  of  load  carrying  capacity. 
Therefore,  postbuckling  analysis  is  essential  to  predict  structural  failure.  We  also  need  a 
failure  criterion  such  as  critical  stress  or  critical  energy  release  rate  to  identify  the  failure 
load. 
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End  shortening  (mm) 


Figure  5.17  Load  versus  end-shortening  with  various  debond  ratios. 


Figure  5.18  Load  versus  out-of-plane  deformation 
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5.6  Comparison  of  Energy  Release  Rate 

First,  in  order  to  verify  the  three  methods  of  computing  G,  using  SEDM,  VCCT, 
and  proposed  CFTM,  a  double  cantilever  beam  modeled  plate  with  offset  node  elements 
was  considered.  The  dimension  of  DCB  specimens  (see  Figure  5.19)  were  that  used  by 
Raju  et  al  [87]  and  are  as  follows:  total  length  101.6mm;  delamination  length  50.8  mm; 
width  24.4  mm;  total  specimen  thickness  3.3  mm;  sub-laminate  thickness  1.65  mm.  The 
material  properties  used  are  listed  in  Table  5.7.  In  addition  to  unidirectional 
graphite/epoxy,  a  16-layer  angle  ply  laminate  with  the  lay-up  [+45,  -45  and  an 
isotropic  DCB  were  also  analyzed.  The  transverse  force  applied  to  each  ligament  of 
DCB  is  1  N/m. 

The  normalized  energy  release  rate  values  computed  using  the  average  G  values 
for  the  various  specimens  are  shown  in  Figures  5.20  to  5.23,  respectively.  The  average 
values  used  for  normalization  were  compared  with  those  of  the  3-D  analysis  results  of 
Raju  et  al.  [87]  in  Table  5.8.  The  average  G  values  from  CTFM  are  closer  to  those  from 
3-D  than  the  average  G  values  from  VCCT.  The  accuracy  of  CTFM  depends  on  the 
refinement  of  the  finite  element  along  the  crack  front  while  the  accuracy  of  VCCT 
depends  not  only  the  refinement  of  the  finite  element  along  the  crack  front  but  also  that 
of  ahead  and  behind  crack-tip.  This  may  be  the  reason  why  the  average  values  of  G  from 
VCCT  deviate  from  those  of  3-D.  However,  the  comparison  of  normalized  G 
distribution  along  the  crack  front  is  very  good  for  the  example  considered. 


108 


Table  5.7  Elastic  material  properties  [87] 


Resin 

Aluminum 

Graphite/Epoxy 

E,  GPa 

3.4 

71.0 

134 

E2  GPa 

3.4 

71.0 

13.0 

Gi2  GPa 

1.3 

27.3 

6.4 

Vl2 

0.3 

0.3 

0.34 

Table  5.8  Average  strain  energy  release  rates  for  DCB  specimens  (10"^  J/m^ ) 


Isotropic 
(aluminum) 

Graphite/Epoxy 
(0  degree) 

Graphite/Epoxy 
(+/-  45  degree) 

Graphite/Epoxy 
(90  degree) 

3-D 

1.0 

0.572 

2.54 

N/A 

CFTM 

0.98 

0.541 

2.71 

5.33 

VCCT 

0.90 

0.522 

2.01 

5.14 

Figure  5.19  Deformed  shape  of  DCB  specimen. 
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Figure  5.20  G  distributions  for  isotropic  DCB  specimens. 


G  distributions  for  G/E  (+/-45)  DCB 
specimens 


Figure  5.21  G  distributions  for  graphite/epoxy  (+/-45  degree)  DCB  specimens. 


110 


G  distributions  for  G/E  (0  deg.)  DCB 
specimens 
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Figure  5.22  G  distributions  for  graphite/epoxy  (0  degree)  DCB  specimens. 


Figure  5.23  G  distributions  for  graphite/epoxy  (90  degree)  DCB  specimens. 
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Next,  the  energy  release  rates  of  stiffened  panel  with  skin-stiffener  debond  from 
Reference  [86]  was  computed  under  the  same  loading  and  boundary  conditions  in 
Reference  [86]  where  stretching  displacements,  25.4x10'^  mm.  were  prescribed  to  the 
model.  Figure  5.24  shows  the  deformed  shape  of  debonded  stiffened  panel  using  3-D 
solid  elements.  On  the  other  hand,  deformed  shape  of  the  debonded  stiffened  panel  using 
shell  elements  is  shown  in  Figure  5.25.  Energy  release  rates  along  the  debond  front  was 
computed  using  VCCT  with  various  mesh  refinement,  and  were  compared  with  those  of 
Reference  [86]  in  Figure  5.26.  The  comparisons  show  that  the  computed  G  distributions 
from  both  the  shell  elements  and  the  3-D  solid  elements  are  similar  to  each  other. 
However,  the  computed  G  from  shell  elements  is  about  20  %  lower  than  that  from  3-D 

solid  elements.  Wang  and  Raju  [85]  proposed  a  procedure  to  reduce  this  discrepancy. 

However  It  is  still  n 

Finally,  under  the  assumption  that  debond  fronts  move  along  the  length-wise 
direction  of  stiffener,  distributions  of  energy  release  rate  along  debond  front  with 
debond  ratios,  a/L=0.35  and  0.5,  were  obtained  from  virtual  crack  closure  techniques  (as 
shown  in  Figures  5.27-5.29).  It  should  be  noted  that  aforementioned  assumption  is  valid 
when  local  buckling  is  dominant  buckling  mode.  As  the  load  increases,  the  energy 
release  rate  underneath  the  stiffener  increases  rapidly.  This  comes  from  the  fact  that  the 
bending  stiffness  of  the  flange  underneath  the  stiffener  is  greater  than  that  of  the  flange 
away  from  the  stiffener.  Thus,  the  differences  of  displacements  and  rotations  for  a  pair 
of  nodes  between  skin  and  flange  underneath  stiffener  is  much  larger  than  between  those 
away  from  the  stiffener. 


112 

Comparison  of  energy  release  rates  versus  end-shortening  using  strain  energy 
derivative  method  and  virtual  crack  closure  technique  was  made  in  Figure  5.30.  In 
general,  energy  release  rate  from  strain  energy  derivative  method  agrees  well  with 
average  energy  release  rate  from  virtual  crack  closure  technique.  The  critical  energy 
release  rate  must  be  obtained  by  experiment.  If  the  assumed  critical  energy  release  rate  is 
about  Gcr  =200  J/m^,  then  debond  extension  will  be  first  started  underneath  the  blade 
stiffener  below  the  buckling  load.  From  the  fact  that  distributions  of  energy  release  rate 
in  the  debond  front  are  not  uniform,  using  the  energy  release  rate  from  virtual  crack 
closure  technique  may  give  more  accurate  estimation  of  debond  extension  than  that  from 
the  strain  energy  derivative  method. 

5.7  Conclusion 

Buckling  and  postbuckling  behavior  of  a  stiffened  panel  with  a  partial  skin- 
stiffener  debond  are  investigated  using  the  finite  element  method.  The  present  model 
shows  good  correlation  with  the  results  of  existing  buckling  analyses  of  a  delaminated 
plate.  There  exists  a  critical  debond  ratio  such  that  if  the  ratio  of  debond  is  less  than  a 
critical  value,  then  the  buckling  load  is  unchanged.  Furthermore,  this  critical  debond 
ratio  depends  on  geometry,  boundary  conditions  and  material  properties. 
Three  methods  for  computing  the  energy  release  rate  along  the  crack  front  with  plate 
elements  are  proposed  to  predict  debond  extension  during  postbuckling.  When  local 
buckling  mode  is  dominant,  the  maximum  energy  release  rate,  which  occurs  below  the 
stiffener  blade,  can  be  much  greater  than  the  average  energy  release  rate. 
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Figure  5.24  Deformed  shape  of  debonded  stiffened  panel  using3-D  solid  elements. 


Figure  5.25  Deformed  shape  of  the  debonded  stiffened  panel  using  shell  elements 
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Through-the-width  position  of  stiff  ener 


Figure  5.27  Distribution  of  Energy  Release  Rate  (a/L=0.34) 
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Figure  5.29  Distribution  of  energy  release  rate  in  width  direction  of  stiffener  (a/L=0.5) 
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Figure  5.30  Comparison  of  energy  release  rate  with  end-shortening  using  strain  energy 
derivative  method  and  virtual  crack  closure  technique. 


CHAPTER  6 
CONCLUSIONS  AND  FUTURE  WORK 


The  major  objective  of  this  study  was  to  analyze  buckling  and  delamination  of 
composite  stiffened  panels  subjected  to  axial  compression. 

First,  analytical  models  using  several  structural  analysis  models  were  used  to 
assess  the  adequacy  of  the  design  model  and  the  correlation  with  experimental  results 
for  a  stiffened  panel  designed  using  the  PASCO  program.  Of  the  effects  neglected  by  the 
simple  model,  shear  deformation  was  the  most  important,  accounting  for  about  11% 
difference  in  buckling  load.  The  effect  of  simplified  (simple  support)  boundary 
conditions  was  small.  The  addition  to  the  design  model  of  shear  loads  and  imperfections 
to  improve  the  robustness  of  the  result  did  help,  even  though  the  inclusion  of  one-sided 
imperfection  apparently  induced  sensitivity  to  imperfection  of  the  opposite  sign.  Overall, 
the  simplified  model  did  produce  designs  that  in  the  experiments  failed  slightly  above 
the  design  load.  Furthermore,  one  of  important  finding  from  the  analytical  and 
experimental  correlation  study  of  stiffened  panel  is  that  an  approximate  modeling  of 
skin-stiffener  intersection  with  common  nodes  can  predict  prebuckling  response  of  the 
stiffened  panel  with  non-negligible  error. 

The   most   significant   difference   between   the   analytical   predictions  and 
experimental  measurements  was  the  substantial  out-of-plane  pre-buckling  deformations. 
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To  explain  these  differences,  imperfections,  load  eccentricities,  and  loading  platen  tilt 
angles  were  considered.  Of  these  the  loading  platen  tilt  produced  similar  patterns  of 
deformation,  but  these  had  more  nonlinear  characteristics  than  the  measured  deformations. 

Next  buckling  and  postbuckling  behavior  of  a  stiffened  panel  with  a  partial  skin- 
stiffener  debond  were  investigated  using  the  finite  element  method.  The  present  model 
shows  good  correlation  with  the  results  of  existing  buckling  analyses  of  a  delaminated 
plate.  There  exists  a  critical  debond  ratio  such  that  if  the  ratio  of  debond  is  less  than  a 
critical  value,  then  the  buckling  load  is  almost  unchanged.  This  critical  debond  ratio 
depends  on  geometry,  boundary  conditions  and  material  properties. 

The  average  G  values  from  CTFM  are  closer  to  those  from  3-D  analysis  than  the 
average  G  values  from  VCCT.  The  accuracy  of  CTFM  depends  on  the  refinement  of  the 
finite  element  along  the  crack  front  while  the  accuracy  of  VCCT  depends  not  only  the 
refinement  of  the  finite  element  along  the  crack  front  but  also  that  of  ahead  and  behind 
crack-tip.  This  may  be  the  reason  why  the  average  values  of  G  from  VCCT  deviate  from 
those  of  3-D  models.  However,  the  agreement  between  the  three  methods  of  normalized  G 
distribution  along  the  crack  front  is  very  good  for  the  example  considered.  . 

When  local  buckling  mode  is  dominant,  the  maximum  energy  release  rate,  which 
occurs  below  the  stiffener  blade,  can  be  much  greater  than  the  average  energy  release  rate. 

Future  work  in  this  area  could  be  in  developing  progressive  damage  models  for 
stiffened  composite  panels  as  well  as  application  of  CTFM  on  delamination  propagation. 
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